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ABSTRACT 


A numerical optimization technique is used to fully automate the trajectory design 
process for an asymmetric configuration of the proposed Advanced Launch system 
(A.L.S.). The objective of the A.L.S. trajectory design process is the maximization of the 

vehicle mass when it reaches the desired orbit. 

The trajectories used in this thesis were based on a simple shape that could be described 
by a small set of parameters. The use of a simple trajectory model can significantly reduce 
the computation time required for trajectory optimization. 

A predictive simulation was developed to determine the on-orbit mass given an initial 
vehicle state, wind information, and a set of trajectory parameters. This simulation utilizes 
an idealized control system to speed computation by increasing the integration time step. 

The conjugate gradient method is used for the numerical optimization of on-orbit mass. 
The method requires only the evaluation of the on-orbit mass function using the predictive 
simulation, and the gradient of the on-orbit mass function with respect to the trajectory 
parameters. The gradient is approximated with finite differencing. 

Prelaunch trajectory designs were carried out using the optimization procedure. For 
the trajectory shape originally used, the procedure proved to be highly sensitive to the 
initial guess of the optimal solution. To rectify this problem, the trajectory shape was 
modified and this change resulted in a procedure that was robust to the initial choice of the 

guess. 

The predictive simulation is used in flight to redesign the trajectory to account for 
trajectory deviations produced by off-nominal conditions -- e.g., stronger than expected 
head winds. For this purpose, only a single trajectory parameter is modified - the value of 
Qa used in the constant aerodynamic loading portion of the trajectory. 
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Chapter One 


INTRODUCTION 

1.1 Background and Problem 

Both manned and unmanned launch vehicles currently require a large amount of 
planning to design trajectories for each mission flown. This leads to high costs and causes 
the launch system to be inflexible to last minute payload or orbital changes. In addition, 
trajectories for these vehicles are designed based on monthly wind averages rather than 
wind data measured the day of launch. If winds differ too widely from those expected for 
the design, the vehicle might reach orbit with insufficient fuel left for required orbit 
maneuvers. 

A new launch vehicle is being developed by NASA to alleviate these problems. Called 
the Advanced Launch System (A.L.S.), the proposed vehicle should decrease the costs 
associated with mission preparation while allowing the vehicle to be robust to unpredicted 
wind variations. The vehicle will have wind velocity data available to it from 
measurements taken a half hour before launch. The goal of the A.L.S. trajectory design 
has been established by NASA as the maximization of the mass of the vehicle once it has 
reached a desired orbit. This mass will be referred to in this thesis as the "on-orbit mass. 
The primary constraint on trajectory design is that the normal aerodynamic loads on the 
vehicle not exceed design limits. 

Several vehicle models of the A.L.S. have been proposed. This thesis focuses on an 
asymmetrical version submitted by General Dynamics in 1988. This vehicle consists of a 
core stage containing the payload attached to a single booster. The booster contains seven 
engines while the core contains only three such engines. The resulting asymmetry in thrust 
leads to trajectories requiring large pitch angles of attack. 
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An earlier study was conducted by Boelitz 1 . The guidance and control concepts he 
developed are included n the full simulation used for this study. Boelitz restricted vehicle 
motion to the pitch plane by nulling both yaw and roll torques. This thesis will be based on 
the same assumption so that comparisons between the two studies can be made. 

A previous thesis in ascent guidance was done by Corvin 2 . Corvin's thesis had the 
same goals as this study but utilized a different vehicle mode, the single stage-to-orbit 
(SSTO) Shuttle n. Corvin developed a simple trajectory shape that could be described by a 
small set of parameters. Boelitz used this trajectory shape and this thesis will also use this 
shape. 

1.2 Method 

To reduce costly prelaunch preparation time and to enhance system flexibility, the 
trajectory design process should be fully automated. Given a specific vehicle mode, 
payload, and wind information acquired shortly before launch, a mission planner should be 
able to use a computer program that determines to within some tolerance on on-orbit mass, 
the set of trajectory parameters which will maximize the on-orbit mass of the vehicle. That 
program should then be able to design a trajectory based on these parameters and save the 
trajectory in the flight computer's memory for the guidance system to command during 
flight. If the flight computers have the computational capability, then this same program 
could be used to redesign the trajectory in flight. This redesign would allow the vehicle 
trajectory to adapt to winds that differ significantly from the prelaunch measurements. 

A computer algorithm has been developed which meets the objectives outlined above. 
The algorithm utilizes a predictive simulation which calculates the vehicle's on-orbit mass 


1 Boelitz, F.W., "Guidance, Steering, Load Relief and Control of an Asymmetric Launch Vehicle". 
1989. Massachusetts Institute of Technology Master of Science Thesis, CSDL Report T-1036. 

2 Corvin, M.A., "Ascent Guidance for a Winged Boost Vehicle." 1988. Massachusetts Institute of 
Technology Master of Science Thesis, CSDL Report T-1002. 
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and a numerical optimization scheme which uses the on-orbit mass to define it s objective 
function. 

The predictive simulation is a simplified simulation that is used in place of the full 
simulation to calculate on-orbit mass. The predictive simulation idealizes the control 
system used in the full simulation. This idealized system greatly reduces the total 
computation time for the on-orbit mass calculation because the vehicle's equations of 
motion can be integrated using a much larger integration time step. The inputs to the 
predictive simulation include the current vehicle stage, the prelaunch wind measurement, 
and a set of trajectory shape parameters. 

The numerical optimization scheme uses the negative of the on-orbit mass as the 
objective function it will minimize. This is equivalent to maximizing on-orbit mass. 
Whenever the optimization scheme is required to evaluate its objective function at a given 
set of trajectory parameters, it utilizes the predictive simulation. 

A conjugate gradient method was chosen for the numerical optimization scheme 
because it requires only the objective function and the gradient of the objective function 
while showing a high speed of convergence. The gradient of the objective function, the 
derivatives of the on-orbit mass function with respect to the trajectory parameters is 
approximated with finite differencing. 

1.3 Overview 

The vehicle model used in this thesis is discussed in Chapter 2. The configuration, 
thrust modelling, mass properties, and aerodynamic characteristics are all presented. The 
environmental modeling, including atmospheric pressure, density, and wind velocity is 
discussed. The chapter also describes the kinematics and dynamic equations of motion 
needed to simulate the vehicle's motion in a pitch plane about a spherical Earth. Finally, 
the chapter describes the normal aerodynamic load constraint placed upon the A.L.S. by the 
vehicle designers. 

Chapter 3 describes the trajectory design, guidance, and control concepts developed by 
Boelitz and implemented in the full simulation used for this thesis. At the end of the 
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chapter, the proposed trajectory optimization procedure is presented and the requirements 
for such a procedure are described. 

The predictive simulation is discussed in Chapter 4. The simplified kinematics and 
dynamics are presented along with the idealized control approximation. The predictive 
simulation was compared to the full simulation for several different time steps and the 
accuracy was very good. 

Chapter 5 justifies the choice of the conjugate gradient method for the numerical 
optimization of on-orbit mass. The underlying theory and procedure is also presented. 

The simulation results are presented in Chapter 6. The numerical optimization 
procedure worked for the trajectory shape initially defined by Corvin but a small 
modification of the shape was made to improve the procedure’s robustness to arbitrary 
initial guesses of the optimal set of trajectory parameters. 

Chapter 7 presents conclusions drawn from the thesis and suggestions for future 
research. 
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Chapter Two 


VEHICLE DESCRIPTION AND MODELING 


2.1 Physical Description of A.L.S. Vehicle 

The vehicle model used for this ascent guidance study is based upon a configuration of 
the Advanced Launch System (A.L.S.) proposed by General Dynamics in 1988. The 
A.L.S. is under development by NASA and the Air Force to provide an unmanned boost 
vehicle capable of delivering large payloads to low-Earth orbit. 

The basic arrangement of vehicle components is illustrated by Figure 2.1. The vehicle 
consists of a core stage and a single booster stage arranged in an asymmetric parallel 
configuration. 

Both the core and booster stages have identical, non-throttleable engines fueled by a 6. 1 
mixture of liquid hydrogen (LH) and liquid oxygen (LOX). Each engine has a vacuum 
thrust level of approximately 612,000 lbs. The booster has seven such engines while the 
core has only three engines. The fuel tanks in each stage are identical. The engine nozzles 
are gimbaUed in both the pitch and yaw directions to provide thrust direction control. Since 
the booster has more engines, its fuel tanks will be depleted before those of the core. 

When this occurs, the booster is separated from the core. 

The upper portion of the core contains the payload bay. The diameter of this section is 
larger than the lower section of the core. The lower section of the core has approximately 
the same dimensions as the booster stage. The inertial measurement unit (IMU) is located 
at the base of the core. 

The booster's engines, servos, and fuel lines are contained in a Booster Recovery 
Module (BRM). This is the only recoverable part of the A.L.S. Separation of the BRM 
from the booster occurs approximately twenty seconds after the booster separates from the 
core. Parachutes are then deployed to return the BRM to Earth and recovery is made at sea. 
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Figure 2.1: A.L.S. Configuration 
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The A.L.S. vehicle uses a total of 10 gas generator engines. Several important features 
of the engines are presented in Table 2.1. The engines are non-throttleable, meaning that 
the thrust level cannot be changed during flight. The mass flow rate and vacuum thrust are 
both assumed to be constant. The thrust magnitude of each engine within the atmosphere 
will therefore vary with atmospheric pressure alone. 


NAME 

SPECIFICATION 

Cycle 

Gas Generator 

Propellants 

LOX/LH 

Throttling Rage 

Fixed 

Propellant Flow Rate 

1 ,427 Lbs/sec 

Vacuum Thrust 

612 KLbs 

Weight 

6,744 Lbs 

Inside Diameter 

88.0 in 

Length 

150 in 


Table 2.1: A.L.S. Engine Characteristics 


The thrust direction is the only available control input to the system. Each engine is 
gimballed in both the pitch and yaw planes so that the deflection angle of the engine nozzles 
can be changed during flight The limit on the gimballing capability of each engine is ± 9 . 
Because of the asymmetry of the vehicle, the engines were installed with a 5* cant angle in 
the pitch plane. This design feature allows the vehicle to have a large gimballing capability 
to prevent limiting of the engine nozzle deflection angle. The rate of change of the nozzle 
deflection is limited to 10*/sec. 
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For this study, it was assumed that the thrust for the set of engines in each stage could 
be represented by a resultant thrust vector such that two thrust vectors, Tb and T c , would 
be controlling the vehicle. This assumption is described by Figure 2.2. As shown, both 
thrust vectors are limited to a deflection of 9* from the installed cant angle of 5*. 


Booster, Core, 



NOTES 

1 ) All 10 engines are installed with a 5° cant. 

2) All 10 engines have the same gimballing capability of 
± 9° from installed cant. 

3) Resultant thrust vector of core, Tc. acts through point A. 

4) Resultant thrust vector of booster, Tb. acts through point B. 


Figure 2.2: Thrust Model 


It was further assumed that both thrust vectors will be deflected by the same angle, 6 . 
This angle is computed by the flight control system so that the vehicle can maintain control 
while steering to the trajectory commanded by the guidance system. 
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2.2 Mass Properties 


To simulate the linear and angular acceleration of the vehicle during flight, knowledge 
is required of the time-varying mass properties: mass (m), eg position, and moment of 
inertia. The mass decreases as fuel is expended. As the mass is decreased, the eg position 
and moment of inertia change. It has been assumed in this thesis that the motion of the 
ALS vehicle is constrained to lie in the pitch plane. Therefore, only the pitch plane 
components of eg location (x cg and z cg ) are required and only the moment of inertia about 

the pitch axis (I yy ) is required. 

A calculation of the dry mass properties is computed before launch and is combined 
with the initial fuel mass properties to give a total vehicle mass property calculation before 
lift-off. The mass properties can then be updated continuously during flight since mass 
flow rate is assumed to be constant . 

The dry model used is the same as that developed by Boelitz, based upon 
recommendations from NASA. The vehicle is separated into various geometric solids and 
shells. All components were assumed to have uniform mass density. 

The fuel tanks in both the core and the booster were modelled as hollow shells with the 
fuel inside modeled as a solid cylinder with time-varying length. The engine modules and 
payload bay were modeled as solid cylinders. The dry mass properties are presented in 
Table 2.2 below: 


Vehicle 

m 

x cg 

z cg 

Iyy 

Component: 

(slugs) 

(ft) 

(ft) 

(slug ft 2 ) 

Core 

10,924 

138.0 

0 


Booster 

5,781 

63.5 

0 


TOTAL: 

16,705 

112.2 

-11.1 



Table 2.2: Dry Mass Properties (Datum at base of core) 
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23 Aerodynamic Characteristics 


In order to compute the aerodynamic forces and moment in the pitch plane, it is 
necessary to know the coefficients of aerodynamic normal and axial force. These 
coefficients are vehicle dependent and are functions of both Mach number and angle of 
attack. NASA provided CSDL with updated aerodynamic data in 1989. Lift, drag, and 
moment coefficients were provided for 0<Mach<8 and for angles of attack between -14* 
and +14*. The vehicle simulations used in this study interpolated between these data points 
to obtain coefficients for the current vehicle state. Linear interpolation was used between 
consecutive Mach numbers and cubic spline interpolation was used between consecutive 
angle of attack values. 

2.4 Environmental Conditions 

In addition to the aerodynamic coefficients described above, knowledge is required of 
the atmospheric model (pressure, density, and speed of sound) and of the winds to 
determine the aerodynamic forces. The vehicle simulations used in this study implement 
the equations for atmospheric pressure, density, and speed of sound given in the 1976 US 
Standard Atmosphere. These parameters are all defined as functions of altitude with a 
range of 0 to 282,000 feet. Above this range, the air density and pressure are assumed to 
be zero and the speed of sound is assumed to be the same as the vacuum speed of sound. 

It is assumed in this thesis that the A.L.S. vehicle will not be capable of sensing winds 
during flight The winds will be measured one half hour before launch using the Jimsphere 
radar-tracked balloon system. It is assumed that this is the only wind information that can 
be used in trajectory design, both before launch and during flight 

For this study, NASA provided a wind profile from Vandenberg AFB. A certain 
percentage of the wind profile is used for trajectory design before launch, and a different 
percentage of the wind profile is used for the in-flight simulation. The goal of this variation 
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is to test the guidance system for robustness in the presence of unexpected wind 
disturbances. 

Because this study is limited to the pitch plane, the winds were assumed to lie in the 
pitch plane and act in a direction that is parallel to the local Earth-relative horizontal. Also, 
the wind profiles were linearized using straight line approximations between a finite 
number of data points. The profile used by this study and its linearized approximation is 
shown in Appendix A. In this profile, the winds dissipated above 66,000 ft. 

2.5 Coordinate Frames and Kinematics 

Three reference frames are used in this study to simulate the motion of the vehicle 
about a spherical Earth. They are defined as. 

(1) Inertial Earth-Centered Reference Frame. (X, Y, Z) 

All equations of motion are referred to this non-rotating reference frame. The 
origin is at the center of the earth. The Z axis points through the North Pole. The 
X axis points through zero longitude and the Y axis completes the right-handed 

set. 

(2) Local Geographic Frame: (un, ue, ug) 

The origin is at the center of gravity of the vehicle. The positive ug axis points 
towards the center of the earth. The u N axis lies on the plane formed by the Z axis 
and u G and points north. The u E axis completes the right-handed set. The wind 
directions and all earth-relative angles are calculated within this reference frame. 
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(3) Body-Fixed Frame: (xb, yB. zb) 

The origin of this frame is fixed to the vehicle's center of gravity and assumes 
there will be no rotation of the vehicle about the xb axis. As stated earlier, this 
thesis constrains the vehicle to unrolled motion in the pitch plane. The xb axis is 
parallel to the centerline of the vehicle and points towards the nose cone. The yB 
axis is in the direction of the cross product of ug and xb- The zb axis completes 
the right-handed set. All forces and torques on the vehicle are computed in this 
frame. 



Figure 2.3: Relationship Between Body and Local Geographic Frames 
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Figure 2.4: Inertial and Body Frame Relationship With Pitch Plane 
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The relationship between the body and local geographic reference frames is shown in 
Figure 2.3. Again, because this study is concerned only with the pitch plane dynamics, the 
heading is determined solely from the initial launch azimuth and the bank angle is set to 
zero. The Earth-relative pitch attitude is the only variable of interest. 

The relationship between the inertial and the body reference frames and the pitch plane 
is illustrated in Figure 2.4. The angle of attack with respect to the Earth-relative velocity, 
a, and the sideslip angle, /?, are also depicted in this figure. 

The rotation of the inertial frame into the body frame was described by an Euler angle 
set, (% 0, 0) representing a sequential azimuth, elevation, and bank transformation. This 
transformation is described in Etkin. The rotation matrix obtained using this transformation 
is: 


cos 0 cos *F cos 0 sin *F -sin 0 


[C] = 


sin 0 sin 0 cos sin 0 sin 0 sin 
- cos 0 sin + cos 0 cos 


sin <Z>cos 0 


cos 0 sin 0 cos cos 0 sin 0 sin cos 0 cos 0 
+ sin 0 sin *¥ - sin 0 cos 


( 2 . 1 ) 


The rotation matrix is used to resolve the components of a vector known in the inertial 
frame into the body-fixed frame and vice-versa. If a vector B is known in the inertial frame 
such that: 


B = X] Uja +3 ’/UyI +Z/UZI = Xb UXB + yB t*YB +ZflUzB (2.2) 


but the components xb, yB* and zb are not known, then these components can be 
determined from the following relation: 


XB ) 

r , IX/ ) 

yB 

=[c] bv 

ZB 1 

\z/ 


(2.3) 
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The rotation matrix, [C], is orthogonal so that [C] 1 = [C] T - Therefore, a vector which is 
known in the body-fixed frame can be expressed in the inertial frame with the relation: 


XI 


ZJ 1 

1 U B 


2.6 Dynamics and Rigid Body Equations 


Figure 2.5 shows the angular rate and moment notation used for the body-fixed frame 
where: 


L = rolling moment 
M = pitching moment 
N = yawing moment 


(Or = rate of roll 
(Op = rate of pitch 
(Oy = rate of yaw 



Figure 2.5: Angular Rates in Body Frame 
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Since this study is limited to the pitch-plane motion, it was assumed that the vehicle 
experiences no roll or yaw torques: 


L = N = 0 (2.5) 

The roll and yaw rates are therefore also zero: 

(0 r = COy = 0 ( 2 . 6 ) 

and the sideslip angle, /?, is zero. 

2.6.1 Flight Orientation Parameters 

The following flight orientation parameters, pictured in Figure 2.6, are used to describe 
the state of the vehicle in the trajectory plane: 

0 = earth-relative pitch attitude 

a = angle of attack with respect to the air-relative velocity 

(Xe = angle of attack with respect to the Earth-relative velocity 

otw = angle of attack contribution from winds = a — gce 

y= flight path angle 

Ve = Earth-relative velocity 

Vw = wind velocity 

Va = air-relative velocity = Ve - Vw 

The air-relative velocity is the difference between the Earth-relative velocity and the 
wind velocity. The angle of attack with respect to the air-relative velocity is equal to the 
sum of the Earth-relative angle of attack and the angle of attack produced by winds. The 
angle of attack contribution from winds is defined as positive if the cross product of Va 
with Ve points in the positive ye direction. The angle of attack with respect to the air- 
relative velocity is used in the calculation of the aerodynamic coefficients. 
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Figure 2.6: Flight Orientation Parameters 
2.6.2 Forces and Torques 

The forces that act on the vehicle during endoatmospheric flight are: the thrust forces 
(Tb and T c ), the aerodynamic forces (F\ and Fa), and the force of gravity (F g ). A free 
body diagram of the vehicle in the pitch plane is given in Figure 2.7. 

The aerodynamic force acts at the vehicle center of pressure and can be resolved into the 
body-fixed coordinate frame. For an unrolled vehicle with zero sideslip angle, the normal 
aerodynamic force acts in the -zb direction and the axial aerodynamic force acts in -xb 
direction. These forces can be written in the body-fixed frame as: 

Fn = - S Q Cm uzb 

F a = -S<2Cauxb (28) 
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where: 


S = reference area = constant 
Q = dynamic pressure = p V* 2 / 2 
p = air density 

Va = magnitude of air-relative velocity 

Qv = coefficient of aerodynamic normal force = Cfj (a, Mach) 
Ca = coefficient of aerodynamic axial force = Ca(oc, Mach) 
uzb = unit vector in z-direction of body-fixed frame 
uxb = unit vector in x-direction of body-fixed frame 



Figure 2.7 : Vehicle Free Body Diagram in the Pitch Plane 
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The aerodynamic pitching moment is expressed in a similar form: 
MaerO = s QCn ilcpx - x cg) + S QCa (-lepz + z cg) 


( 2 . 9 ) 


where: 


l cpx = location of center of pressure with respect to datum 

/ cpz = location of center of pressure with respect to centerline of core 

For this study, both l cpx and l cpz were assumed to be constants. 

The two thrust forces, T b and T c , are expressed in the body-fixed frame as: 

T b = T b cos 8 uxb + Tb sin 8 uzb 

T c = T c cos 8 uxb + T c sin 8 uzb 


where: 


Tb = booster thrust 
T c = core thrust 
8 = nozzle deflection 

The thrust contribution to the pitching moment exerted on the vehicle is: 

Mthrust = T b x cg sin 8 + T c x cg sin 8 -T b (D + z cg ) cos 8 - T c z cg cos 8 ( 2 . 12) 


where: 

x cg = body x-axis eg position (measured from datum at base of core) £ 0 
z cg — body z-axis eg position (measured from centerline of core) ^ 0 
D = distance between centerlines of core and booster = constant 
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The force of gravity points toward the center of the Earth along uq and can be 
expressed in the inertial frame as: 


Fg = - m g (R //?) 


(2.13) 


where: 


R = position vector of vehicle expressed in inertial coordinates 
R = magnitude of inertial position vector 
m = m(r) = vehicle mass at time t 


2.6.3 Rigid Body Equations of Motion 

The forces acting on the vehicle can be summed in the body-fixed frame to give the 
resultant force in body-fixed coordinates: 

Fnet = Fn + Fa + Tb + T c + Fg (2.14) 

This force can be resolved into the inertial frame by the use of equation 2.4. The 
translational equations of motion are then: 

<LK = y 

dt (2.15) 


<UL-(-L) F 

dt m net (2.16) 


The rotational equations of motion are calculated using the angular momentum 
principle: 


M 


dt relative to 

body frame 


0) x [/]© 


(2.17) 
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(2.18) 


This set of equations can be expanded into the following: 

M = ^co + [/}^ + ©x[/]to 

where M is the total moment acting on the vehicle and is expressed in the body-fixed frame 
as: 

M = {m) (2.19) 

and co is the angular velocity of the body expressed in the body-fixed frame as. 

- M 

CO = [(Op) 

\cOyl (2-20) 


and [/] = [/(r)] = vehicle inertia at time t 

To simplify equation 2.18, two assumptions can be made. First, the body-fixed axes 
are assumed to form a principal axes set for the vehicle so that the inertia matrix, [/], is 
diagonal. Second, the term that consists of the product of the time derivative of inertia with 
angular rate was found to form a very small contribution to the total moment and was 
therefore neglected. Using these assumptions, equation 2.18 can be written in scalar form 

as: 

L = IxX^-0>pCOy{lyy-Izz) ( 2 . 21 ) 

M = Iyy^-(Oy(Dr{I Z i-Ixx) (2.22) 

N = /„ ^ - COrCOp {I xx - Iyy) (2.23) 
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These equations are known as Euler’s equations of motion. However, we stated earlier 
that: the rolling moment, yawing moment, rate of roll, and rate of yaw are all equal to zero: 

L =N - a>r = CDy =0 (2.24) 


Therefore, equations 2.21 and 2.23 are both trivial and equation 2.22 simplifies to: 



(2.25) 


where the pitching moment, M, acting on the vehicle is the sum of the contributions from 
the thrust forces and from the aerodynamic forces: 

M = M A ero + fd thrust (2.26) 


The Euler angle rates are then calculated from the body rates by the use of the following 
non-orthogonal rate transformation matrix: 


[W] = 


1 

0 

0 


sin <t> tan 0 
cos 0 
sin 0 sec 0 


cos d>tan 6 
-sin 0 

cos 0 sec 0 . 


(2.27) 


Using this transformation matrix, the time rate of change of the Euler angle set which 
describes the attitude of the vehicle with respect to the inertial frame is given by: 


<L 

dt 


0 

0 


= [W] 


(Or 
(O p 
a )y 


(2.28) 


where: 


(Of = C0y = 0 for all time t (assumption made for pitch plane analysis) 
top = body pitchrate = dd/dt 
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In summary, the aquations of motion implemented in the six degree of freedom 
simulation used for this thesis are. 

Translational: 

fLR = V (2.29) 

d t 

(LV=(X)F (2.30) 

di m NET 
Rotational: 

lilt = 7,J M (2.31) 

dx 7y 


<*> 


'(Or' 

0 

= [W] 

(Op 

.T'- 


Uy. 


(2.32) 


2.7 Constraints 

The normal aerodynamic force can produce a large bending moment on the vehicle.as 
the vehicle moves through the atmosphere. Structurally, them is a limit on the magnitude 
of the bending moment that the vehicle can sustain without failure. Therefore, it ts 

necessary to constrain the normal aerodynamic force during flight 

Recalling equation 2.7, the magnitude of the normal aerodynamic force can be 

expressed as: 
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where: 


S = reference area = constant 
Q = dynamic pressure = p V^ 2 / 2 
p = air density 

V* = magnitude of air-relative velocity 

Cs - coefficient of normal aerodynamic force = Cn (at, Mach) 

For small values of angle of attack, the coefficient of normal aerodynamic force can be 
approximated as a linear function of angle of attack. Given this approximation, the 
coefficient can be expressed as: 

Cn = Cno « ( 2 . 34 ) 


where: 


Cno — 


&Cn 

da 


The normal aerodynamic force can then be approximated by: 

Fn *= S Q Csa Q 


( 2 . 35 ) 


( 2 . 36 ) 


Since S is a constant and C^ a is approximately equal to a constant, the normal 
aerodynamic force is roughly proportional to the product of dynamic pressure, Q, and the 
angle of attack, a. Therefore, to control the normal aerodynamic force, the vehicle is 
usually controlled to limit the product of Q and a. A limit on Qa has been specified by the 
designers of the A.L.S. vehicle. This limit is the primary constraint on the vehicle's 
trajectory during endoatmospheric flight. 
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The dynamic pressure is proportional to the product of air density, p, and the square of 
the magnitude of the air-relative velocity, V A . The air-relative velocity monotonically 
increases from zero during flight while the air density monotonically decreases with 
altitude. The result of these effects is that the dynamic pressure rises from zero to a 
maximum within the atmosphere and then decreases back to zero when the vehicle has left 
the Earth's atmosphere. A typical dynamic pressure profile showing this behavior is 
illustrated in Figure 2.8. 

There are several methods which have been developed to control the Qa product so that 
the constraint on normal aerodynamic force is not exceeded. Corvin, in his Shuttle II 
study, used the throttling capability of that vehicle's engines to control velocity and thus 
dynamic pressure. Boelitz, in the initial A.L.S. study, used an angle of attack limiting 
mode within the control system. Both Corvin and Boelitz found that appropriate trajectory 
shaping could also help to keep the Qcc product below the specified limit. 



Figure 2.8: Typical Dynamic Pressure Profile 
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Chapter Three 


TRAJECTORY DESIGN, GUIDANCE, AND CONTROL 

CONCEPTS 


3,1 Introduction 

The new trajectory design and guidance concepts developed in this study are based 
upon those developed by Boelitz. Trajectory design is the process of choosing a trajectory, 
based upon objectives and constraints, for the guidance system to command during flight. 
Boelitz’s objective for the trajectory design of the guidance commands was to achieve as 
high an "on-orbit mass" as possible while constraining the vehicle to fly within the designer 
specified Qa limit. The "on-orbit mass" is defined as the mass of the vehicle when it has 
just been inserted into the desired, elliptical orbit. An advantage of maximizing on-orbit 
mass is that it allows more fuel to be available for any post boost maneuvers. Also, 
minimizing the fuel needed for endoatmospheric boost will allow the A.L.S. to carry larger 

payloads into orbit. 

A disadvantage of the Boelitz study was that the guidance was implemented in an open- 
loop manner such that the trajectory was not updated during flight Consequently, if wind 
dispersions from the pre-launch measurement occur or if an engine-out failure occurs, then 
the trajectory designed before flight might not be optimal for on-orbit mass. In this thesis, 
an attempt is made to update the trajectory during flight using the vehicle's current state, a 
predictive simulation for on-orbit mass, and a numerical optimization scheme to maximize 
for on-orbit mass. The control and estimation concepts used in this thesis are the same as 
those used in the previous study. 

This chapter reviews the work done by Boelitz and introduces the improvements that 
are made in trajectory design and guidance. Section 3.2 describes the mission of the 
vehicle and how its flight is divided into separate phases. Each phase necessitates the use 
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of different guidance and control schemes. Section 3.3 describes which signals can be 
sensed and which signals must be estimated. This section then reviews the estimators that 
Boelitz used for angular velocity, angle of attack, dynamic pressure, and acceleration 
direction. The guidance and control concepts used in flight are discussed in Section 3.4. 
Section 3.5 reviews the pre-launch trajectory design profiles and methods for each phase. 
Finally, Section 3.6 introduces the method for automating the pre-launch trajectory design 
and updating the trajectory in flight. 

3.2 Mission and Flight Phases 

The A.L.S. must be able to carry payloads through the atmosphere to a low-Earth 
parking orbit. The amount of lead-time and planning for the mission must be kept to a 
minimum. The vehicle must also have the capability to tolerate wind dispersions from pre- 
launch measurements and still have enough fuel left after atmospheric ascent to maneuver 
into the desired parking orbit. 

The ascent trajectory of the A.L.S. can be divided into four distinct phases. Each phase 
involves different constraints and thus different guidance and control schemes. A 
schematic of the ascent is shown in Figure 3.1. 

Phase 1 is a vertical rise from the launch pad so that the vehicle can clear the launch 
tower. The pitch attitude of the vehicle during this phase must be held constant at 90°. 
Attitude control is thus used to meet this objective. 

Phase 2 is characterized by a rapid pitchover designed to orient the vehicle to the initial 
state required by Phase 3. The state required at the start of Phase 3 is very different from 
the state achieved at the end of Phase 1 . Therefore, Phase 2 must quickly pitch over the 
vehicle while the dynamic pressure, Q, is still small enough so that the Qoc product remains 
well below the designer specified Qa limit. The guidance profile used in this phase is 
based on pitch attitude rate calculated as an analytical function of time. 
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At=7 


Figure 3.1: A.L.S. Flight Phases 
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Phase 3 is the endoatmospheric phase during which the trajectory may be constrained 
by the Qa limit. The Qa product reaches its maximum during this stage and could nse 
above the specified Qa limit if it were not constrained. During this phase, the guidance 
command profile is a stored acceleration-direction trajectory. A dual-mode control system 
is used that will follow the acceleration-direction profile if the Qa product is below the 
specified limit, but will switch to a Qa limiting control mode if the Qa product is above the 

limit. 

Phase 4 begins when the atmospheric density has fallen to a sufficiently small value so 
that the dynamic pressure is low once again. The transition from endoatmospheric to 
exoatmospheric flight occurs during this phase. A modification of the Space Shuttle 
Powered Explicit Guidance system (PEG) is used to guide the vehicle. Because Phase 4 
starts in the upper atmosphere, the PEG algorithm neglects aerodynamic forces. At the end 
of this phase, the vehicle is inserted into low Earth orbit. The specified parameters for this 

orbit are: 


(1) Radius of perigee = r pe = 80 nautical miles 

(2) Radius of apogee = r ap = 150 nautical miles 

(3) Horizontal velocity at perigee: 

1/2 


f perigee 


2f* I, 

Vpe 


L r op + r P e J 


where p = gravitational constant 


3 3 In-Flight Guidance and Control 

This section presents a review of the guidance and control concepts developed by 
Boelitz in the previous A.L.S. study. This thesis involves modifications of the guidance 
but does not attempt to improve on the control concepts. The control loops and estimators 
used in this study are therefore the same as used in the previous study. 
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33.1 Phase 1 and 2 Guidance and Control 

Figure 3.2 shows a general block diagram describing the guidance and control for both 
Phase 1 (vertical rise) and Phase 2 (launch maneuver). The guidance commands of both 
Phase 1 and Phase 2 are based on pitch attitude rate, co c . For Phase 1 the vehicle must 

rise vertically to clear the tower and so: 


o>c(t)= 0 

<W= 90 ’ 

0c(O = 0c lNrr = 90- 


(3.1) 


For Phase 2, the vehicle must start at the final conditions for Phase 1 and then orient 
itself to meet the starting conditions for Phase 3. Corvin and Boelitz both used a sinusoidal 
pitch attitude rate for this purpose whose form is given by the following relation: 


Ofc(r) = Q | 1 - COS (t- T Vert) 


Tvert ^ ^ Tzick + Tv„ 


(3.2) 


where: 


Q = half the maximum pitch rate 

Tvert = duration of Phase 1 (time required for vehicle to clear launch tower) 
T Kick = duration of Phase 2 


The shape of this maneuver is shown in Figure 3.3. Both f2 and TKick are constants 
determined by the trajectory design process discussed in the next section. With this 
analytical function, the integration of 0) c shown in Figure 3.2 can be performed analytically 

giving the following relation for commanded pitch attitude: 


6 c (t) = (t- Tvert) - ^2^* s ’ n k ^ v ' erf )]) + ^ ClNtT 


Tv.rtSliTtid+Tv.r, ( 3 . 3 ) 
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Referring again to Figure 3.2, the guidance employed in both Phase 1 and Phase 2 is 
open-loop, meaning that the guidance commands to the control system are designed prior to 
launch and are not updated during flight. This study did not attempt to close the guidance 
loop during Phase 1 and 2 because it was decided that the altitude of the vehicle during both 
of these phases was sufficiendy low so that wind dispersions were not a serious problem. 

The attitude commands generated by the open-loop guidance system are fed directly 
into the control system which tries to null the error between the commanded attitude and the 

sensed attitude, e«. The attitude error is compensated by proportional-integral control. The 

resulting signal is then summed with the fed-forward commanded pitch rate, at, to produce 
a net pitch rate command, of. This signal is compared to the estimated pitch rate resuldng 

in a pitch rate error signal, 

This error signal is multiplied by a proportional gain, Kt/Kv, to provide an engine 
nozzle angle command to the engine nozzle servos. Boelitz linearized the vehicle 
dynamics, using the assumption of a planar trajectory, to find a transfer function between 
engine nozzle deflection, & and pitch attitude, ft On the basis of this linearization, the gain 
of this vehicle dynamics transfer function is approximately Kv where: 

T x C o 

Kv = -j^ (3.4) 

lyy 
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This is a quantity which monotonically increases with time. Therefore, Ky, is used in the 
denominator of the proportional gain to keep the total inner loop forward gain constant 

For this thesis, the engine nozzle servos were idealized such that the commanded nozzle 
deflection was perfectly achieved: 

8=8 C (3.5) 

33.2 Phase 3 Guidance and Control 

Traditionally, launch vehicles have used "acceleration direction" steering for this phase 
of flight. The "acceleration direction" is the inertial acceleration of the vehicle excluding the 
acceleration caused by gravity. The control loop is commanded to follow either an 
acceleration direction or a pitch attitude command from the guidance system. In addition to 
the nominal acceleration direction (or attitude) and attitude rate feedback signals in the 
control loop, a parallel "load relief signal is added. This add-on load relief employs a 
signal proportional to the measured AV component acting normal to the vehicle’s 
longitudinal axis. At low frequencies, this signal is approximately proportional to the angle 
of attack. The effect of the add-on load relief is to rotate the vehicle in the direction of the 
air-relative velocity vector and thereby reduce the angle of attack. The disadvantage of the 
add-on load relief concept is that the load relief is in conflict with the acceleration direction 
control and can produce significant deviations from the desired trajectory in the presence of 
winds. The control concept used in this thesis is one which overcomes some of the 
disadvantages of the traditional load-relief concept. It is a dual mode, acceleration 
direction/(2a-limiting guidance and control scheme that was first developed by Bushnell. 
Boelitz later applied the concept to the A.L.S. vehicle in his study. 

A general block diagram for Phase 3 acceleration direction/Qa-limiting guidance and 
control is shown in Figure 3.4. The control system has two modes. In the primary mode, 
the vehicle's "sensed" acceleration is commanded to follow a stored acceleration direction 
profile. The sensed acceleration includes only the contributions to vehicle acceleration from 
thrust and aerodynamic forces. Gravity forces are not included. The secondary mode (the 
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Qa limiting mode) is only activated when the vehicle experiences large aerodynamic loads. 
These loads are constantly monitored by comparing the predicted angle of attack, a pre d, to 
an angle of attack limit defined by the division of the Qa limit by the estimated dynamic 

pressure, Q: 

Opred ~ ® + e A (3.6) 


Glim — 


I QG\lim 

/V 

Q 


(3.7) 


where: 

a = estimated angle of attack 
e A = acceleration-direction error 

If the predicted angle of attack is greater than the angle of attack limit, then the control 
system will switch to the secondary mode and thus null out any difference between the two 
signals. Using this dual-mode concept, acceleration-direction trajectory following is 
unimpaired by an add-on load relief function and the vehicle autopilot will only perform 
load relief, in the form of Qa -limiting, when needed. 

The option for either closed-loop or open-loop guidance is shown on the block diagram 
in Figure 3.4. For closed-loop guidance, a new acceleration direction profile can be 
provided by trajectory design in flight and a new Qa -limit can be used for the angle of 
attack limit calculation. This closed-loop guidance concept will be discussed further in 
section 3.6. 
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Open-Loop 

Guidance 



Figure 3.4: Phase 3 Guidance and Control Block Diagram 
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The control system shown to the right of the mode switch in Figure 3.4 is similar in 
form to that used for Phase 1 and 2. The switch provides either an acceleration direction 
error signal, e A , or an angle of attack error, e a , to the control system. The error signal 
used will be modified by proportional-integral control producing a commanded pitch 
attitude rate, co c '. This signal is compared to the estimated pitch rate, 0). The resulting 
error in pitch rate, is multiplied by a time- varying proportional gain, Kg/Ky, to provide 
an engine nozzle deflection command, 8 C , to the engine nozzle servos. Again, the nozzle 
servos are idealized such that the commanded nozzle deflection is perfectly achieved. 

Because the two control modes have different dynamic characteristics, different sets of 
control gains must be used. The alternative sets of gains are designated by the two 
different subscripts in the gains K Pi l , K hx and Kg l 2 shown in Figure 3.4. Boelitz, using 

a linear approximations for the system, showed how the two modes involved different 
dynamics and did a stability analysis for each. He also reset the integrator error so that the 
nozzle command would be continuous when switching between modes. 

3.4 Sensing and Estimation 

3.4.1 Sensed Signals 

It is assumed for this thesis that the Inertial Measurement Unit (IMU) contains an 
attitude measuring device that measures vehicle attitude with respect to an inertial reference 
frame. The measured attitude is defined by a set of Euler angles. In addition, the IMU 
contains an integrating accelerometer that measures "sensed" inertial velocity. This velocity 
is the integral of sensed acceleration- i.e. acceleration produced by thrust and aerodynamic 
forces and excluding gravity. The IMU signals are processed to yield the following signals 
required for guidance and control: 

1) Vehicle pitch attitude relative to a local earth horizontal: 6 

2) Inertial velocity increments: AV 
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3) Inertial velocity: V 

- this signal is the sum of the sensed inertial velocity and the integral of gravity, 
where gravity is specified by a gravity model 

4) Inertial position: R 

- the integral of inertial velocity 

In addition to these four processed signals, the nozzle deflection, 5, is determined from 
nozzle actuator measurements. 

There are several signals needed for the guidance and the control of the vehicle that 
cannot be measured directly during flight. These are the pitch angular rate (co), the angle of 
attack with respect to the air-relative velocity (a), the dynamic pressure ( Q ), and the 
acceleration direction of the vehicle excluding gravity (ffo). All of these signals must be 
estimated during flight. 

3.4.2 Angular Rate 

The pitch angular rate estimate, co, is used as the inner loop feedback variable during 
flight phases 1, 2, and 3. It is also used for the estimation of angle of attack, a to account 
for the effects of tangential and centripetal acceleration at the IMU. Because there can be 
significant quantization and intrinsic noise in the IMU attitude measurement, it is not 
sufficient to use a derived rate signal based only on the quotient of pitch attitude change and 
the sampling period. The noise in the attitude signal will be magnified in the derived 
attitude rate, especially if the sampling period is small. 

Boeltiz used an angular rate estimator based on a first order digital complementary filter 
which has both derived rate (change in pitch attitude over a sampling interval) and estimated 
angular acceleration as inputs. A continuous-time representation is shown in Figure 3.5. 
Derived rate is the low frequency input and is passed through a low-pass filter to produce a 
rate estimate that is accurate at low frequencies. In the implementation of the high 
frequency path, an equivalent representation is used in which an estimate of angular 
acceleration is passed through a low-pass filter. This procedure gives the same results as 
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high-pass filtering a high frequency rate signal based upon integration of angular 
acceleration. The angular acceleration estimate is based upon the IMU-measured inertial 
velocity, the engine deflections, and corrections for any small errors in the modelling of the 
thrust and aerodynamic forces. 


low- pass filter 


Otow frequency 


(£) ► CO 


G>hlgh frequency 


TS 


xs + 1 



1 


xs + 1 



high-pass filter 


Figure 3.5: Continuous Signal Representation of Angular Rate Estimator 
3.43 Angle of Attack 

The angle of attack estimate, a, is employed in Phase 3. It is used in the mode 
switching logic and it is the primary outer-loop feedback variable in the control system of 
the Qa limiting mode. Boeltiz used a second order digital complementary filter for this 
estimator. The continuous-time representation of the filter is shown in Figure 3.6. The 
low frequency input, ai ow frequency, is an angle of attack estimate determined from the 

following procedure: 

1) Estimate the normal acceleration at the center of gravity from the IMU sensed 
velocity increments and estimated angular rate and angular acceleration. 

2) Subtract the thrust contribution to the estimated normal acceleration to isolate the 
component produced by normal aerodynamic force. 

3) Estimate the magnitude of the normal aerodynamic force by multiplying the estimated 

normal acceleration produced by this force with the vehicle mass. 
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4) Determine the normal aerodynamic force coefficient by dividing the estimated normal 

aerodynamic force by the product of the reference area and the estimated dynamic 
pressure. 

5) Use the normal aerodynamic force coefficient and Mach number to search through 

the aerodynamic data tables to find a corresponding value of angle of attack. 


low-pass filter 



Figure 3.6: Continuous Signal Representation of Angle of Attack Estimator 

The high frequency input to the complementary filter is based upon pitch attitude. 
Referring back to Figure 3.6, it can be seen that an s multiplying the high frequency input, 
othigh frequency, is equivalent to using dhigh frequency as the high frequency input into a lower 
order filter. The equivalent filter is shown in Figure 3.7. 

The problem now lies in estimating a discrete-time othigh frequency- As shown in 
Chapter 2, the angle of attack can be expressed as: 

a=6-y+a w (3.8) 
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low-pass filter 



Figure 3.7: Alternate Representation of Angle of Attack Estimator 


Taking the derivative of both sides yields: 


d=d-y+ aw 


(3.9) 


This equation can be expressed in the discrete-time domain as: 

Aa = A6 - Ay+ Aaw (3.10) 

The incremental change Aaw is cause by variations in the winds normal to the velocity 
vector. Because these variations cannot be measured, this quantity is neglected in the angle 
of attack estimation. The quantity Ay is also neglected to yield the relationship: 

A« = A6 0.11) 

The omission of the incremental change in flight path angle can cause a small time- 
varying bias in the estimation of Aa. A second order complementary filter is used for the 
angle of attack estimation so that this bias can be attenuated. 
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3.4.4 Dynamic Pressure 


Phase 3 utilizes estimated dynamic pressure, Q, in both the mode switching logic and 
in the Qa limiting mode. The dynamic pressure estimate is also used in the angle of attack 
estimator. 

The dynamic pressure is a function of both atmospheric density and the magnitude of 
the air-relative velocity of the vehicle. The atmospheric density is assumed to be known 
during flight based upon a standard atmospheric model. However, because the vehicle has 
no wind sensors, the magnitude of the air-relative velocity must be estimated. Boelitz used 
the earth-relative velocity measurement from the IMU along with the estimated angle of 
attack and the assumption of horizontal winds to accomplish this task. 

3.4.5 Acceleration Direction 

The estimated acceleration direction, 0 A , is used as a feedback signal in the acceleration 
direction mode of Phase 3. The only measurements that are used to estimate this quantity 
are the inertial velocity increments processed from the IMU accelerometer measurements. 
Each control cycle, the inertial velocity increments are expressed in the body frame as: 

AV^ = increment in velocity along the vehicle x (roll) axis 

AV 2 = increment in velocity along the vehicle y (pitch) axis. 

AV 3 = increment in velocity along the vehicle z (yaw) axis. 

The direction of the acceleration vector is then computed in terms of the pitch and yaw 
angles, (3 p and P y : 

P p = tan 1 (-4V 3 MV,) (3.12) 

P y = tan* 1 (AV 2 IAV x ) (3.13) 
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The AV measurements are noisy signals because of IMU quantization. Therefore, the pitch 
and yaw acceleration angles are each sent through a low-pass filter with the following 
continuous-time form: 


Pis) _ 1 

P(s) r P S+l ( 3 - 14 ) 


where x p is the filter time constant The unit vector, U A > representing the estimated filtered 
acceleration direction in body axes is calculated as. 


Ua = Unit value of 



(3.15) 


The estimated acceleration direction angle in the pitch plane is then: 


6a = tan 1 



(3.16) 


In the simulation of the acceleration direction control system, the error signal between 
the commanded acceleration direction, d Ac , and the estimated acceleration direction is 
determined by the following procedure. The cross product between Uac (the commanded 
acceleration direction unit vector) and U a is calculated. 

C = U A x Uac (3.17) 


The angle, p A , between the two vectors is: 

)3 A = sin -1 1 C | (3.18) 

The vector composed of the error angles in roll, pitch, and yaw is then. 

U E = [ unit (C) ] (3.19) 
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The pitch component of this vector is the error signal, e a, shown in Figure 3.4. 

3.5 Pre-Launch Trajectory Design 

3.5.1 Introduction 

A guidance system which is implemented in an open-loop manner in flight will utilize 
commands that were determined prior to launch by a trajectory design scheme. In the 
previous A.L.S. study, the guidance for phases 1 through 3 was open-loop. The guidance 
for Phase 4 was calculated closed-loop by a version of the Powered Explicit Guidance 
(PEG) program used in the current Space Shuttle system. This program neglects the 
effects of aerodynamic forces and can only be used after the vehicle has reached the upper 
atmosphere. 

Suitable trajectories must meet the desired objectives while satisfying the specified 
constraints. For the A.L.S. vehicle, the primary objective of trajectory design is to reach 
orbit with as much fuel left as possible. The vehicle should have a high on-orbit mass, 
mfinai . A large on-orbit mass will allow for any post-boost maneuvers that might be 
needed. Minimizing the fuel needed for endoatmospheric boost would also allow for 
heavier payloads to be carried into orbit. The primary constraint on the trajectory design is 
that the Qa product experienced during flight must be within the designer specified limit at 
all times. 

The trajectory design process is greatly simplified if a trajectory shape or form is 
specified for each phase. The trajectory shapes for the first three phases are described by: 

Phase 1: The altitude the vehicle must reach at the end of the vertical rise. 

Phase 2: The sinusoidal function of pitch attitude rate versus time. 

Phase 3: The three parameters of an angle of attack profile. 

The remainder of this section is devoted to the explanation of the trajectory shapes and the 
pre-launch trajectory design concepts that Boelitz developed in his thesis. For the present 
study, the entire pre-launch trajectory design procedure has been automated using a 
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numerical optimization technique. The automation of the pre-launch trajectory design 
process is discussed in Section 3.6. Section 3.7 discusses how this same technique is 
applied to in-flight trajectory design. 

3.5 2 Phase 1: Vertical Rise 

The only parameter that can be varied for this phase is the altitude the vehicle is required 
to reach at the end of the vertical rise. Although this parameter has the possibility of 
affecting the rest of the trajectory, it was decided to hold this parameter constant at 400 feet 


3.53 Phase 2: Launch Maneuver 


This phase is very important to the overall trajectory of the vehicle. It must be able to 
take the vehicle from the final state of Phase 1 to the initial state required for Phase 3. The 
launch maneuver must accomplish this quickly while the vehicle remains well below the 
Qa limit. 

As discussed earlier in this chapter, the trajectory shape used for this phase is based on 
a sinusoidal function of pitch angular rate: 


(Dc = Q[ 1 - COS (f- TverA ) 

\ \T Kick -11 


Tver! - 1 - Tgidt + Tvtn 


( 3 . 20 ) 


where: 


Q = half the maximum commanded pitch rate 
Tvert = duration of Phase 1 
T K^k - duration of Phase 2 


Integration of this equation yields: 


where: 


Qf= final pitch attitude of the vehicle at the end of the launch maneuver 
6i = initial pitch attitude of the vehicle at the start of the launch maneuver = 90' 

It should be noted that because phase 1 involves a vertical rise from the launch tower, 
6i will always be 90*. With this knowledge, the only two parameters that can be adjusted 
to change the shape of the pitch rate profile are 8/ and T^ick- 

A constraint on the choice of these parameters is that the angle of attack reached by the 
end of the launch maneuver, a/, must match the angle of attack desired for the beginning of 
Phase 3. With this constraint in mind, Boelitz developed a procedure that automated the 
design of this phase of the trajectory. He chose 8f and a j as the two inputs to the 
procedure. The procedure used a reduced order/idealized control simulation to predict the 
vehicle state at the end of the launch maneuver. Using this simplified simulation of the 
vehicle's motion, he iterated on the value of Txick until the magnitude of the difference 
between the desired final angle of attack, a/, and the predicted final angle of attack was less 
than a small tolerance. The golden section search in one dimension was used for this task. 

3 .5.4 Phase 3: Angle of Attack Profile 

In Phase 3, the trajectory design process is based upon a simple angle of attack profile. 
As the trajectory design program simulates the vehicle following this angle of attack profile, 
the acceleration direction of the vehicle is calculated and stored away for use in the guidance 
loop during flight. 

The shape of this angle of attack profile is shown in Figure 3.8. Figure 3.9 shows the 
angle of attack guidance and control system used for the trajectory design simulation of this 
phase. The error in angle of attack is sent into the same control system that was used for 
attitude control in phases 1 and 2. The angle of attack is first commanded to be a constant 
value, a \ . This value of angle of attack must be equal to the final angle of attack, a y, 
achieved at the end of Phase 2. The angle of attack limit is continuously monitored by 
dividing a specified Qa limit by the current Q. When this limit becomes less than <X\, the 
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vehicle is commanded to fly along the limit. Because the specified Qa limit is a constant, 
the angle of attack limit is inversely proportional to the dynamic pressure. This results in a 
"bucket" shape for the angle of attack limit. The vehicle can move off of the limit once the 
limit becomes larger than a second constant angle of attack, a 2 . There are therefore three 
parameters which define the shape of this trajectory: a u a 2 , and the specified Qa limit. 
This profile is constrained in several ways: 

1) The Qa limit used must be less than the designer specified Qa limit by a finite 

amount so that the normal loads on the vehicle will not be larger than structural 
limits. 

2) Phase 2 must be able to meet a\. 

3) The need for a 2 to be compatible with Phase 4 requirements so that a "smooth 

transition to Phase 4 can occur. 



Figure 3.8: Phase 3 Angle of Attack Profile for Trajectory Design 
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Figure 3.9: Phase 3 Control System For Trajectory Desgin 
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3.5.5 Phase 4: Powered Explicit Guidance 

Phase 4 is initiated when the vehicle is in the upper atmosphere and aerodynamic forces 
arc at a minimum. For this phase, a version of the Powered Explicit Guidance (PEG) 
program used on the Space Shuttle is used to guide the vehicle into the desired orbit. This 
guidance program is closed-loop in that it produces its own guidance commands in flight. 

It is not necessary to do any trajectory design prior to launch for this phase of flight. 
However, this program is very important to the overall trajectory design because it is used 
to predict the on-orbit mass of the vehicle. The primary objective of the overall trajectory 
design is to produce a high on-orbit mass. 

The method updates the commanded acceleration direction angle in pitch every six 
seconds using a "linear tangent guidance law" of the form: 

tan 6a = Kq + (f -to) (3.22) 

where 6 A is the commanded acceleration direction angle in the pitch plane and Kq, to, and 
K\ are parameters which the program adjusts to minimize the propellant required to bnng 

the vehicle into orbit, thus maximizing on-orbit mass. 

This program treats the vehicle as a point mass so that rotational dynamics are neglected 
and the model is reduced from six degrees of freedom to three degrees of freedom. The 
only variable inputs to this program are the vehicle's inertial position, velocity, acceleration 
and acceleration of gravity vectors. The program is then used to predict the on-orbit mass 

of the vehicle. 

3.6 Automation of Trajectory Design 

The trajectory of the vehicle has been segmented into distinct phases with each phase 
defined in terms of a simple "shape". Using these shapes, the endoatmospheric boost 
portion of flight, phases 1 through 3, can be completely described with the following set of 

trajectory parameters: 
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1) 6f = final pitch attitude of the vehicle at the end of Phase 2 

2) a\ = final angle of attack for phase 2 = initial angle of attack for Phase 3 

3) as = final angle of attack for Phase 3 

4) Qa limit = specified Qa limit for mission < vehicle designer specified Qa limit 

This is a very simple shape because only four parameters are needed to describe the 
trajectory. The problem is then how to choose these four parameters so that the on-orbit 
mass is maximized while the vehicle flies within all of the constraints considered in Section 
3.5. 

Boelitz tuned each of these parameters by trial and error using the full six degree of 
freedom simulation until it seemed that on-orbit mass had reached an maximum. This 
method is time-consuming and does not guarantee that an maximum has been reached. 

In this thesis, a method has been developed which will "automatically" determine the 
set of trajectory parameters which maximize on-orbit mass. This method is based upon a 
numerical optimization scheme which uses on-orbit mass as the objective function it seeks 
to maximize. A predictive simulation is used to calculate the on-orbit mass. A simplified 
schematic of the automated trajectory design process is shown in Figure 3.10. 

The predictive simulation is initialized with the current state of the vehicle. The only 
disturbance information available to the simulation is the pre-launch wind measurement. 
The optimization algorithm supplies the set of trajectory parameters. The predictive 
simulation will integrate the equations of motion from the current time to the time when 
PEG is used to predict on-orbit mass. The optimization algorithm will continue a 
multivariable search for the trajectory parameters which optimize on-orbit mass until a 
maximum on-orbit mass has been found. 
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This scheme is used before launch to reduce launch preparation time and to allow for 
unexpected changes in payload weight. However, because this scheme is "automatic" it 
can be used in flight to update the trajectory in the event of wind dispersions from the pre- 
launch measurement. The predictive simulation stores the complete state history so that a 
new acceleration direction profile can be used for guidance. In flight, there are no 
aerodynamic sensors so the disturbance information is assumed to be the same as it was 
prior to launch. However, if there arc wind dispersions, these will affect the current state 
of the vehicle. This updated state information can then be used to determine a new set of 
flight parameters for the calculation of an updated acceleration direction profile. 

The following two chapters describe the components of this scheme in detail. Chapter 
4 discusses the predictive simulation. A reduced order/ idealized control simulation is used 
to speed computation by increasing the integration time step. Chapter 5 discusses the 
numerical optimization scheme chosen for this problem. A conjugate gradient method is 
used because the problem is highly nonlinear. The gradient is approximated with finite 
differencing. 
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Chapter Four 


PREDICTIVE SIMULATION 


4.1 Introduction 

In Chapter 3, a procedure was described for selecting trajectory parameters to maximize 
on-orbit mass for specified values of the Qa limit. In order to automate this process, a 
numerical optimization scheme is employed for which the objective function is the on-orbit 
mass. The prediction of on-orbit mass is computationally expensive because the vehicle's 
equations of motion must be integrated from the given initial conditions to the desired end 
conditions. A simplified predictive simulation has been written which greatly reduces the 
amount of computation needed to calculate the on-orbit mass. This chapter describes the 
characteristics of this simulation and compares its performance with the full six degree of 
freedom simulation (which includes the control systems described in Chapter 3). 

Section 4.2 justifies and details the reduced-order model used for the predictive 
simulation. In order to increase the integration time step of the simulation, an idealized 
control system was developed. This approximation is described in Section 4.3. Section 
4.4 describes the program flow and presents results for various wind conditions. Remarks 
and conclusions are given in Section 4.5. 

4 2 Reduced-Order Model 

The vehicle's pitch dynamics can be reduced to only three degrees of freedom for this 
study since the vehicle motion is restricted to the trajectory pitch plane. By choosing the 
inertial y axis such that it is perpendicular to the trajectory plane, only one angle is needed 
to specify the orientation of the vehicle. In addition, only two translational variables are 
needed to describe the location of the vehicle center of gravity with respect to the origin of 
the inertial frame. For the predictive simulation, it is still necessary to calculate mass 


- 54 - 



properties, aerodynamic data, environmental conditions, and thrust in the same way as for 
the 6DOF simulation, but the coordinate frame kinematics and equations of motion of the 
vehicle are greatly simplified by choosing the inertial frame in this manner. 

The two reference frames used in the predictive simulation are shown in Figure 4.1 and 
are defined as: 

(1) Inertial Reference Frame: (x,y, z) 

The equations of motion used in the predictive simulation are referred to this 
reference frame. The origin is fixed to the surface of the flat Earth at the launch site. 
The z axis points toward the center of the Earth and the x axis points downrange. The 
y axis completes the right-handed set. 

(2) Body-Fixed Frame: (xb, yB, zb) 

The origin of this frame is fixed to the vehicle's center of gravity and assumes no 
roll motion. The xb axis is parallel to the centerline of the vehicle and points towards 
the nose cone. Because the vehicle is restricted to pitch rotation only, the ya axis is in 
the same direction as the inertial y axis. The zb axis completes the right-handed set. 

The motion of the body-fixed frame with respect to the inertial frame is constrained to 
translation in the inertial x-z plane and rotation about the inertial y axis. The position 
vector, R, locates the origin of the body-fixed frame with respect to the inertial origin. The 
rotation of the inertial frame into the body frame is described by a single angle 8 , the pitch 
attitude. The rotation matrix that transforms a vector from inertial to body coordinates is 
given by: 

cos 8 0 -sin 8 

[^predictive sim = 0 1 0 

. sin 8 0 cos 8 J (4.1) 
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Figure 4. 1 : Predictive Simulation Coordinate Frames 

The aerodynamic and thrust forces are defined the same way in the predictive 
simulation as they were in the full 6DOF simulation: 

F n = - S Q C N u ZB (4 2) 

¥ a = -SQC a u X b (43) 

Tb = Tb cos 8 uxb + Tb sin 8 uzb ^ 4 ' 4 ^ 


T c = T c cos 8 uxb + T c sin 8 uzb 


(4.5) 
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The angle that the vehicle subtends around the Earth during boost is very small and so it 
is possible to approximate the Earth as being flat during boost. Using this approximation, 
the gravity vector will always point in the direction of the positive inertial z axis: 

F g = m g u z (4.6) 

The aerodynamic and thrust forces are resolved from the body-fixed frame into the 
inertial frame by using the inverse of the rotation matrix defined in equation 4.2. These 
forces are then summed with the gravity force in the inertial frame to give the net force 
acting on the vehicle in inertial coordinates: 

FNET = F N + F A +T b + T c + F g (4.7) 

Since it is assumed that there are no out-plane-forces acting on the vehicle, the y- 
component of the net force is always equal to zero: 

Fnet = Fnet x ux + (0)uy + Fnet z uz (4.8) 


The translational equations of motion have the same form as for the 6DOF simulation: 


t£R = v 

d t 


(4.9) 


dX. 

dt 


= (J-)F 

m NET 


(4.10) 


where V = inertial velocity = Ve 

Because there are no out-of-plane forces, the y-component of both position and velocity is 
always equal to zero: 


R = R x u x +(0) uy + Rzuz 

(4.11) 

V = V x u x + (0) uy + Vz uz 

(4.12) 
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The total moment acting on the vehicle is the same as was derived in Chapter 2: 

M - (Tb + T c )x cg sin 8 - Tb(D + z cg ) cos 8 - T c z cg cos 8 

+ S Q Cn ( l C px - x C g) + S Q Ca {'lepz + z cg) (4. 13) 

The rate of pitch is then derived from: 

djpL = / - 1 Af (4.14) 

dt 

The pitch attitude, 0, is the only angle needed to specify the attitude of the vehicle and is 
defined by: 

=(O p (4.15) 

di F 

4 3 Idealized Control 

Boelitz used the full 6DOF simulation for pre-launch trajectory design and tuned the 
flight parameters by trial and error. When the full 6DOF simulation was used for trajectory 
design, pitch attitude control was utilized for phases 1 and 2 and angle of attack control was 
used for phase 3. The pitch attitude controller and angle of attack controller (described in 
Chapter 3) were digitally implemented in the 6DOF simulation with a sampling time of 0.1 

seconds. 

A larger integration time step is desired for the predictive simulation so that it will have 
a faster computational speed than the 6DOF simulation. To achieve this goal, an idealized 
control system was developed for the predictive simulation that replaces both the pitch 
attitude and angle of attack control systems. The idealized system is designed to 
approximate the low frequency response of the actual control systems. The idealized 
control system has several benefits that significantly reduce computation time. The 
replacement of the control systems with idealized control allows a very large integration 
time step to be used for the simulation. Also, the computations required to determine the 
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actuator commands for the actual control systems are replaced by much simpler 
computations for the idealized system. Finally, since the 6DOF control system always 
commands either pitch attitude or angle of attack the use of idealized control allows pitch 
attitude to either be specified or calculated directly from angle of attack. Thus, the 
rotational equation of motion need not be integrated twice to solve for pitch attitude. 

The predictive simulation is used for trajectory design purposes so it must mimic the 
control systems used by the full 6DOF simulation for trajectory design. For the full 
simulation, the trajectory design procedure uses pitch attitude control for phases 1 and 2. 
Angle of attack control is used for phase 3 trajectory design. If the control systems are 
idealized then the assumption is made that the commanded control variable will be perfectly 
achieved: 


Phase 1: 

e=e c = 90' (4.i6) 

Phase 2: 

e= e c = oi [ (t- Tver,) - Sin (t- Tver,)]) + <Ur (4-17) 

Phase 3: 

/< M 

a= ac= i <Xii m ) (4.18) 

\ «2 I 


Using the idealized control assumption, at least one flight orientation parameter has 
been specified for each phase. However, both pitch attitude and angle of attack must be 
known at all times during flight. Pitch attitude is needed to specify the kinematics of the 
body frame with respect to the inertial frame and the angle of attack is needed to determine 
the aerodynamic coefficients and thus the forces acting on the vehicle. 

Using the relationships between the flight orientation parameters shown in Figure 4.2, 
it is possible to calculate the angle of attack from the idealized pitch attitude and vice versa. 
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Figure 4.2: Flight Orientation Parameters 

The following equation can be derived from Figure 4.2: 

6 = a + 7 - aw 


(4.19) 


where: 


6 = earth-relative pitch attitude 

a = angle of attack with respect to the air-relative velocity (Va) 
cqy = angle of attack contribution from winds — cc— (Xe 
gce = angle of attack with respect to the Earth-relative velocity 
y= flight path angle 
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Both the flight path angle, y, and the angle of attack contribution from winds, aw, can 
be determined from the assumed wind model and the translational equations of motion. 
Only one flight orientation parameter is specified by the idealized control assumption at a 
time, either pitch attitude or angle of attack. The unkown flight parameter can always be 
determined from equation 4.20 or equation 4.21: 

Phase 1 & 2: 0 = 0 C therefore a = 0 - y + Ow (4.20) 

Phase 3: a = a c therefore 6 = a + 7 - ccw (4.21) 

The 6 DOF control systems of phases 1, 2, and 3 were used to determine the nozzle 
deflection, S, necessary to rotate the vehicle to meet either the commanded attitude or angle 
of attack. With these control systems idealized, the question becomes how to deflect the 
engine nozzles to meet the assumptions given in equations 4.17 and 4.18. 

It is possible to solve for the nozzle deflection angle from the pitch moment given by 
equation 4.14. This equation can be rewritten using the small angle approximation where 
sin(5) = 8 and cos (8) = 1 : 

M = I6 = (T b + T c ) x cg 8 -T b D - (T b + T c ) z cg + SQC N (l cpx - x cg ) 

+ SQCa {~lcpz **■ %cg) (4.22) 

The above equation can now be solved for nozzle deflection: 

_ _ T b D •+ T Zcg - SQCj y ( lcpx ~ ~ 8QCa (• lepz Zqg) + 10 

T x cg (4.23) 

The pitch moment, 10, is the only term in the above equation that cannot be computed 
in flight from knowledge of the current vehicle state. The vehicle's pitch moment is known 
analytically for the first and second phases of flight using the idealized pitch attitude control 
assumption: 

Phase 1: 10 = I0 C = 0 ( 4 - 2 4) 
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(4.25) 


Phase 2: 16 = I6 C = 




The vehicle's pitch moment during the third phase has a very small average value 
because the pitchrate is only slowly varying after the launch maneuver. Therefore, for the 
purpose of determining nozzle deflection, it was decided to approximate the pitch moment 
of the vehicle as zero during phase 3: 

Phase 3: 16 = 0 (4>26) 


4.4 Predictive Simulation Flow and Results 

A flowchart of the predictive simulation is shown in Figure 4.3. The predictive 
simulation is initialized with the current vehicle state. The vehicle can be either on the 
launch pad or somewhere in flight. The wind information given to the predictive 
simulation will always be the wind profile measured prior to launch. The vehicle and 
environmental calculations are the same as those used in the 6DOF simulation. These 
computed variables include mass properties, atmospheric density and pressure, 
aerodynamic coefficients, winds, and kinematic transformations. 

The guidance procedure commands pitch attitude during phases 1 and 2 and angle of 
attack during phase 3. The idealized control routine then assumes perfect control of one of 
the two flight parameters, attitude or angle of attack, and computes the unknown flight 
parameter. The pitch moment is specified (zero for phases 1 and 3, nonzero for phase 2) 
and is used to calculate the nozzle deflection needed to achieve perfect control. The net 
acceleration is then calculated in the inertial reference frame and the equations of motion are 
integrated using the fourth-order Runge-Kutta method. All of the vehicle, environmental, 
guidance and control calculations are updated in the middle of the time step. The simulation 
continues until the Powered Explicit Guidance (PEG) routine is called upon to predict the 
on-orbit mass, m f , at 120 seconds. At this time, the aerodynamic forces are small in 

comparison to the thrust forces. 
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Figure 4.3: Predictive Simulation Flow Chart 
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To facilitate the comparison between the two simulations, it was decided to use an 
equatorial trajectory in the 6DOF simulation. This eliminated a complicated kinematic 
transformation of the predictive simulation's position, velocity, and acceleration into an 
inclined trajectory plane which would have been required for comparison of results. The 
relationship between the inertial reference frames of the 3DOF predictive simulation and the 
6DOF simulation is shown in Figure 4.4. Because both frames are inertial, the 
transformation between the two frames is constant. The position, velocity, and acceleration 
vectors of the predictive simulation are expressed in the inertial coordinates of the 6DOF 
simulation at the transition to PEG. 



Figure 4.4: Relationship Between Inertial Reference Frames of Full (6DOF) and 

Predictive (3DOF) Simulations 
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The objective of the predictive simulation is to give an accurate prediction of the 
vehicle's flight while keeping computation to a minimum. To illustrate the performance of 
the predictive simulation, the following states were compared to those produced by the full 
6 DOF simulation: 

(1) angle of attack (a) 

(2) flight path angle (y) 

(3) height (H) 

(4) nozzle deflection angle (5) 

Two sets of comparisons were run. The first set corresponded to starting both the 
6DOF simulation and the predictive 3DOF simulation at t = 0 (i.e. at launch) and running 
both until t = 120 seconds. This set will be referred to as the "entire" boost set because the 
predictive simulation flies the entire trajectory. The second set of comparisons was made 
by initializing the predictive simulation with the vehicle state given by the 6DOF simulation 
at t = 60 seconds. The predictive simulation then flew until t = 120 seconds. Thus, the 
predictive simulation flew only a "partial" boost trajectory. This test was done to 
demonstrate improvement in accuracy of the predictive simulation when it is initialized with 
the vehicle state at a later time in flight. 

A study was made to select a value of integration time step for the predictive simulation. 
It was observed that the performance was greatly enhanced by using a small time step, dt, 
of 0.1 seconds for phases 1 and 2. This is due to the fact that the dynamics of the vehicle 
during the launch maneuver are much faster than during the rest of the flight. For phase 3, 
four different time steps were considered: 

dt = {0.1, 0.3, 0.5, 1.0} seconds (4.27) 

The first set of comparison runs corresponding to "entire" boost is shown in Figures 
4.5 to 4.8. The entire state history is shown. For each plot, the solid line represents the 
state variable computed by the 6DOF simulation. The four dashed lines represent the state 
variable computed by the 3DOF simulation for the four different time steps given above. 
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Gamma (deg) Alpha (deg) 



Figure 4.5: Angle of Attack Comparison Between Full and Predictive Simulations 

For Entire Boost 



Figure 4.6: Flight Path Angle Comparison Between Full and Predictive Simulations 

For Entire Boost 
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Nozzle deflection (deg) Height (feet) 



Figure 4.7: Height Comparison Between Full and Predictive Simulations 

For Entire Boost 



Figure 4.8: Nozzle Deflection Comparison Between Full and Predictive Simulations 

For Entire Boost 
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The arrow on the plots represents the direction of increasing integration time step for 
the predictive simulations. 

Figure 4.5 shows the comparison in angle of attack between the simulations. The 
correspondence between the predictive simulation and the 6DOF simulation during phases 
1 and 2 was very good. During phase 3, the 6DOF simulation utilizes proportional-integral 
control for angle of attack while the 3DOF utilizes idealized control. Even though the 
control is idealized for the predictive simulations, errors occur during the Qa limiting 
because the angle of attack limit is computed as the quotient of the Qa limit and Q. The 
dynamic pressure, Q , is a function of the air-relative velocity and this variable is calculated 
from the integration of the translational equations of motion. Increasing the time step of the 
predictive simulation causes integration errors to build in the air-relative velocity. The 
errors in air-relative velocity, through dynamic pressure, can thus enter into the calculation 
of the angle of attack limit. The predictive simulation will always steer perfectly to the 
angle of attack limit, but this limit may not be the same as for the 6DOF simulation because 
of the error in air-relative velocity. 

It can also be seen in Figure 4.5 that the angle of attack given by the 6DOF simulation 
falls in between the angle of attack predicted by the 3DOF simulations for time steps of dr = 
0.3 and 0.5 seconds. Thus, the angle of attack predicted by the 3DOF simulation with dt = 
0.5 seconds is more accurate than the angle of attack predicted by the 3DOF simulation 
with dt = 0.1 seconds. This is a result of the fact that there are two main sources of error 
resulting from the use of the predictive 3DOF simulation. These error sources tend to 
produce angle of attack errors of opposite sign. The idealized control produces a negative 
error in angle of attack and the use of large time steps produces a positive error. For time 
steps between .3 and .5 seconds the errors tend to be offsetting. 

Figures 4.6 and 4.7 illustrate the errors in flight path angle and height. The height for 
the 6DOF simulation also falls within the height curves generated by the 3DOF simulations. 
The errors between the 3DOF simulations and the 6DOF simulation build steadily with time 
because of integration errors in the equations of motion. 
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Figure 4.8 illustrates the error in engine nozzle deflection. The nozzle deflection shown 
includes the initial cant of 5*. The nozzle deflection of the 6DOF simulation has distinct 
spikes whereas the nozzle deflections of all of the 3DOF simulation runs are missing these 
spikes. The spikes in the 6DOF occur around the times when the wind profile has a 
discontinuity in slope. The vehicle must compensate for these discontinuities by rapidly 
pitching over. The pitch moment during these times is thus larger than normal. As a 
result, the zero pitch moment approximation used for phase 3 of the idealized control in the 
3DOF simulation is not valid at these times. The average error of the nozzle deflection in 
the 3DOF simulations grows with time yet is still very small. 

Table 4.1 shows the absolute errors in the final states of the predictive simulations at 
time =120 seconds. The absolute error in predicted on-orbit mass is presented along with 
the percentage error of fuel left in the core on-orbit. This quantity is calculated as follows: 


% fuel/ error = 

} f Ue lftDOF 


(4.28) 


where: 

fuelffiooF = core fuel on-orbit for 6DOF sim = m^DOF ~ m DRY 
m DRY ~ dry mass of core 


3DOF dt 

(sec) 

a error 

(deg) 

yerror 

(deg) 

//error 

(feet) 

terror 

(deg) 

nif error 
(slugs) 

% fuelf 
error 

0.1 

0.61 

1.55 

2,850 

0.16 

7.0 

5.0 

0.3 

0.25 

0.83 

1,390 

0.07 

2.8 

2.0 

0.5 

0.06 

0.15 

20 

0.01 

6.9 

mm 

1.0 

0.60 

1.43 

3,260 

0.17 

15.3 

10.9 


Table 4.1: End State Error Comparison Between Full (6DOF) and Predictive (3DOF) 

Simulations After Entire Boost 




































The second set of comparison runs corresponding to "partial" boost is shown in 
Figures 4.9 to 4.12. As stated above, this set of runs was made by initializing the 
predictive simulation at the state given by the 6DOF simulation at time = 60 seconds. The 
accuracy of the predictive simulations is improved because the modeling and control errors 
have less time to grow. This improvement in accuracy is demonstrated in Figure 4.11 
where the difference between simulations cannot even be observed over the time history of 
height. The absolute errors in end states for this set of runs and the on-orbit mass errors 
are shown in Table 4.2. The errors in the end states for the partial boost runs are 
significantly smaller than those that were produced in "entire" boost runs. Also, the 
percentage error of the core fuel left on-orbit is improved for most of the time steps. 



Figure 4.9: Angle of Attack Comparison Between Full and Predictive Simulations For 

Partial Boost 
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Figure 4.10: Flight Path Angle Comparison Between Full and Predictive Simulations 

For Partial Boost 



Figure 4. 1 1 : Height Comparison Between Full and Predictive Simulations 

For Partial Boost 
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Nozzle Deflection (deg) 



Figure 4.12: Nozzle Deflection Comparison Between Full and Predictive Simulations 

For Partial Boost 


3DOF dt 

(sec) 

a error 

(deg) 

yerror 

(deg) 

//error 

(feet) 

terror 

(deg) 

n* error 

(slugs) 

% fuelf 
error 

0.1 

0.09 

1.48 

630 

0.0277 

3.1 

2.2 


0.03 



0.0035 



0.5 

0.04 



0.0002 

SB 

3.9 


0.17 

1.11 


0.0328 


5.8 


Table 4.2: End State Error Comparison Between Full (6DOF) and Predictive (3DOF) 

Simulations After Partial Boost 
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4.5 Conclusions 


The predictive (3DOF) simulation developed for this thesis eliminates the kinematics 
and dynamics for the more complete simulation. The idealized control assumption also 
significantly reduces computation time and allows the use of a larger integration time step. 
These simplifications and assumptions greatly reduce the computational load of the 
predictive simulation. The accuracy of the predictive simulation has been shown to be very 
good in comparison to the 6DOF simulation. 
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Chapter Five 


NUMERICAL OPTIMIZATION 


5.1 Introduction 

The primary objective of trajectory design for the A.L.S. is the maximization of on- 
orbit mass. An optimization procedure is described in this chapter which will maximize on- 
orbit mass. This mass is determined by the trajectory that is flown and is, therefore, a 
function of the particular parameters chosen to specify the trajectory shape. For a given set 
of trajectory parameters, the on-orbit mass is determined by using the predictive simulation 
described in the previous chapter. The simulation approach is required because it is not 
practical to develop a closed-form solution for the on-orbit mass. Consequently, the on- 
orbit mass function is not analytical. Any optimization procedure which seeks to maximize 
on-orbit mass must, therefore, be numerical in form. For the approach used in this thesis, 
the on-orbit mass function is multidimensional because as many as three parameters are 

used to define the trajectory shape. 

This chapter first compares in Section 5.2 several multi-dimensional numerical 
optimization algorithms that are described in the current literature. The particular method 
chosen was a version of the conjugate gradient method. The overall procedure for 
implementing the conjugate method and some of the underlying theory is described in 
Section 5.3. Sections 5.4 and 5.5 describe separate subroutines that had to be performed 
in conjunction with the algorithm. Section 5.4 describes how the on-orbit mass function 
was optimized along a specific search direction. Section 5.5 describes how the gradient of 
the on-orbit mass function was approximated using finite differencing. 
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5.2 Comparison of Numerical Optimization Methods 


There are many multi-dimensional numerical optimization algorithms that have been 
developed over the years. To determine the method most suitable for the optimization 
problem of this thesis, the current literature was studied. Two primary sources were used: 
Scales 3 and Press et. al. 4 . 

All numerical optimization schemes rely on function information and sometimes 
derivative information. To understand these methods and the differences between them, it 
is important to provide some background information and notation relating to the definition 
of a multi-dimensional function and its derivative information. An objective function, F(x), 
which a function of rt independent variables must be defined. These variables form the 
column vector x of length n. All of the optimization methods studied were designed to 
minimize the objective function. The objective of A.L.S. trajectory design is the 
maximization of on-orbit mass. Since the maximization of on-orbit mass is equivalent to 
the minimization of the negative of on-orbit mass, the objective function will be defined as: 

F(x) = - ntf (5.1) 

where m/ is the on-orbit mass of the A.L.S. and is determined from the predictive 
simulation. There are three parameters that define the shape of the trajectory in the 
predictive simulation (Of, C(\ t and a 2 ) and so the vector x is: 

x = [6/ a, cc 2 ]T (5.2) 


3 Scales, L.E., Introduction to Non-Linear Optimization . 1985. London: MacMillan Education 


LTD., pp. 1-106. 

4 Press, W.H., Flannery, B.P., Teukolsky, S.A., Vetterling, W.T., Numerical Recipes . 1986. 
Cambrdige: Cambridge University Press., pp. 274-311. 
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The gradient vector of an objective function, g, is comprised of the n first partial 
derivatives of the objective function evaluated at x: 

g = VF(\) (5 ' 3) 

and provides information on the shape of the function. For the objective function based 
upon on-orbit mass, this gradient is not available analytically and must be approximated 
with finite differencing. This procedure is discussed in Section 5.5. 

Further information on the shape of a function is provided by a matrix composed of the 
n 2 second partial derivatives of F(x) which is called the Hessian matrix. 

G(x) = V 2 F(x) < 5 - 4 ) 


and can be represented by the following tensor notation: 



d 2 F 

dxidxj 


(5.5) 


This matrix is also not analytical if on-orbit mass is used to define the objective function. 
To approximate this matrix with finite differencing would require a large number of 
function evaluations and would be computationally expensive. 

The objective function is used in all numerical optimization schemes. The gradient 
vector and Hessian matrix can also be used in the optimization process because they 
provide information about the shape of the objective function. It was found that the current 
numerical optimization schemes can be divided into three main groups based upon the kind 
of function information they utilize. 

The members of the first group are commonly called "direct search methods . These 
methods use only function evaluations. No information about the gradient vector or 
Hessian matrix is utilized. Without any knowledge of the shape of the function, these 
methods may take an excessively long time to converge to a minimum. 
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An improvement to the direct search methods is provided by a second group of 
algorithms that utilize first derivative information in addition to function evaluations. 
Specifically, they require knowledge of the gradient vector, g. Included in this category are 
the method of steepest descent, the conjugate gradient methods and variable metric 
methods. Using the derivative information, these algorithms are able to change the vector x 
in a direction that will tend to minimize the objective function. The first derivative 
information makes these methods more efficient than direct search methods because these 
methods tend to avoid examining any points that would increase function value. Direct 
search methods are useful for highly discontinuous functions where the gradient vector can 
be singular. However, for the trajectory problem in this thesis, the on-orbit mass is not a 
discontinuous function. For continuous functions, with continuous first derivatives, the 
gradient methods show faster convergence than direct search methods. 

The third group requires information about the Hessian matrix of the objective function 
in addition to the gradient vector. Thus, more information about the shape of the function 
is utilized. This group includes Newton's method and variations of Newton's method. 
These methods show fast convergence for functions with a known Hessian matrix, but the 
cost in computer time of approximating a non-analytical Hessian matrix and the 
inaccuracies that could be introduced through such an approximation makes this method 
impractical for the optimization of on-orbit mass. 

A member of the second group has been chosen for this study because it is possible to 
approximate the gradient of the on-orbit mass function using finite differencing. Of this 
group, the method of steepest descent was eliminated from consideration because its 
convergence is usually slower than the conjugate gradient method. The comparison 
between these two methods is presented in the following section. Both the conjugate 
gradient method and the variable metric methods require only gradient information. The 
variable metric methods construct a rough approximation to the Hessian matrix with 
gradient information collected over a successive number of iterations. Of these two 
methods, the conjugate gradient method was chosen over the variable metric methods for 
its simplicity. 
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53 Conjugate Gradient Method 

Gradient methods, like other search methods, are all iterative: they begin with an initial 
guess of the vector x that will produce a minimum of the objective function, and iterate on x 
until the function is minimized to within some tolerance. The iteration on x is achieved by 
the following relationship: 

xjt+i = x* + 77* P* 

where x* is the previous estimate of the state which will minimize the objective function 
and x*+i is the updated estimate. The vector p* is defined as the search vector , and Tj k is 
a positive scalar weighting which is chosen to give the greatest decrease in the function 
F(x), along the search vector. The quantity, r \ k , is found by a separate one-dimensional 

optimization procedure described in Section 5.3. 

In the method of steepest descent, the search vector is defined as: 


where g* is the gradient vector. The steepest descent method is useful for a function 
where the gradient vector always points in the direction of the minimum. Such a function 
is shown in Figure 5.1 where the function contours are circular in the region around the 
minimum. In this case, the method of steepest descent would find the minimum in one 
iteration since all gradient vectors point toward the minimum. In general, this does not 

occur. 

Using the steepest descent method, a new point, x*+y, is computed when an r\ k is 
found which minimizes F(x) along p*. The gradient at this point, VF(x*+y), is then 
calculated and is orthogonal to p*. The new search vector, p*+y, will therefore also be 
orthogonal to the previous search vector, p*. If the function contours happen to be 
elliptical, as illustrated in Figure 5.2, then the method of steepest descent leads in a slowly 
converging zig-zag pattern to the minimum. 
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Figure 5. 1 Steepest Descent Path for Circular Function Contours 

The conjugate gradient method makes an improvement on the definition of the search 
vector direction and thus shows faster convergence than the method of steepest descent in 
most cases. To understand the properties of this method, it is useful to examine the case 
where the objective function to be minimized is a quadratic function of n independent 
variables: 

F(x) = x r Ax + b 7 x + c (5.8) 

where A is a constant symmetric n by n matrix, b is a constant n-vector and c is a scalar. 
The gradient of this function is: 

g(x) = Ax + b (5.9) 

The Hessian matrix for the quadratic function is simply: 

G(x) = A = a constant matrix (5.10) 
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Figure 5.2 Steepest Descent Path for Elliptical Function Contours 

If A, and thus G, is positive definite, then the quadratic function will be convex and have a 
global minimum with no local minima or saddle points. For the rest of the discussion, it is 
assumed that the quadratic function defined in equation 5.8 does have a positive definite A 
matrix. 

A unique property of quadratic functions can be developed by taking the gradient at a 
point xi: 

g(xi) = Axi+b < 511 > 

and again at a point X 2 : 

g(x 2 ) = Ax 2 + b < 5 - 12 > 
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Subtracting equation 5.1 1 from equation 5.10 yields: 


g(xi) - g(x 2 ) = Axi - Ax 2 = A (xj - x 2 ) (5.13) 

Substituting equation 5.10 into the above equation results in the following property which 
is unique to quadratic functions: 

g(xi)-g(x 2 ) = G(x]-x 2 ) (5.14) 

Having defined the above property, it now possible to describe the conjugacy property 
of the conjugate gradient method for a quadratic function. The conjugate gradient method 
will minimize a quadratic function with a positive definite Hessian matrix in n iterations or 
less. This is accomplished by choosing the search directions used in equation 5.6 such that 
they are mutually conjugate to the Hessian matrix: 

p/Gp* =0 for j *k (5.15) 

This conjugacy condition can be changed into a more physically meaningful form by 
first multiplying both sides of the equation by the scalar Tty: 

rjj p/ Gp* = 0 (5.16) 

Using equation 5.6, this is equivalent to: 

(xy +1 -xy)7"Gp*=0 (5.17) 

Substituting the transpose of equation 5.13 results in: 

(&+i-«/) T P*=0 (5-18) 

The above equation means that the current search direction, p*, must be orthogonal to 
all of the changes in gradients, (gy+i - gy) for j = 0(l)Jfc-l, that were defined previously. 
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This property assures that the minimum of a quadratic function with a positive definite 
Hessian matrix will be found in at most n iterations. A proof of this is given in Scales. 


Scales also derives the following definition of a search vector which satisfies the conjugacy 
condition: 

o (5.19) 

P*= -g k + PkPk-l 


where the initial search vector is: 


There are several versions of the conjugate gradient method and they differ only in the 
definition of /)*• The particular version chosen for this thesis is the Polak-Ribiere 

method where /J*is given by the expression: 

p k = (g k - gfc-i) r g* (5.20) 

g/t-i T g*-i 

The Polak-Ribiere algorithm was recommended by both Scales and Press for its efficiency 
and stability in finding the minimum of a non-quadratic function. 

It should be noted that the conjugate gradient method's defimuon of the search vector 
does not require any more function evaluations or gradient approximations than the steepest 
descent method. The previous gradient vectors are simply stored and used in the calculation 
of the new search vector. 

The conjugate gradient method will only calculate search vectors that are mutually 
conjugate and will converge to a minimum in n iterations if the objective function is 
quadratic with a positive definite Hessian matrix. Most objective functions used in 
numerical optimization procedures, including the on-orbit mass function used in this thesis, 
are not quadratic functions with a known positive definite Hessian matrix. Saddle points 
and local minima could occur. Therefore, the gradient algorithm will generally go through 
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more than n iterations for convergence. If local minima are present in the objective 
function, then convergence to the local minimum is not guaranteed. 

The algorithm flow is shown in Figure 5.3. The routine is supplied with an initial 
guess of the minimum, x 0 . The routine then computes the function and gradient, using the 
predictive simulation, at this initial point. The first conjugate gradient direction chosen is in 
the direction of the negative gradient because no previous gradient information is available. 
The routine then begins the main iteration loop. At the start of each iteration, a separate 
search algorithm (line minimization) finds the scalar, rjt-y, that minimizes VF(x^.y) along 
the search vector direction, Pk-i • Equation 5.6 is then used to update the estimate to x*. 
The method used for this search algorithm is described in section 5.3. The gradient is then 
calculated at this new point and the information is used to establish a new search direction, 
p*. A convergence check is made between the line minimization block and the gradient 
calculation block. The convergence is based on the difference of the function values 
between separate conjugate directions. If this difference is smaller than a specified 
tolerance, then the algorithm returns with the most recent estimate of the minimum. 

5.4 Minimization Along Search Direction 

The one-dimensional minimization of the objective function that must be carried out 
over Tjk for each k-iteration may require many iterations. Recalling equation 5.6: 

X*+1 = X*+ ThcPk < 5 - 22 ) 

The scalar 77 * must be chosen so that the value of the function at x*+i is minimized 
along the search direction p*. This is accomplished as follows: the one-dimensional 
minimization algorithm computes a value for the scalar tj* and a separate routine takes this 
value and substitutes it into equation 5.22 along with the current estimate, x*, and search 
direction, p*. This defines a new estimate x*+i. The function of this new estimate is then 
evaluated and returned to the one-dimensional minimization algorithm. 
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Figure 5.3 Polak-Ribiere Conjugate Gradient Algorithm for Function 
c Minimization 
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The one-dimensional minimization algorithm used in this thesis was comprised of two 
separate programs that are described in Press. The first program brackets the minimum 
with the triplet of abscissas: A, B, and C such that: 

A<B<C (5.23) 


and: 


F(B ) < F(A) 


(5.24) 


and: 


F(B) < F(C) (5.25) 

Equations 5.24 and 5.25 state that, of the three points which bracket a function, B has 
the lowest function value. If the bracketing conditions shown above are achieved, then a 
minimum of the function must exist somewhere within the interval [A, C]. 

Examples that illustrate whether or not a set of three points meet the bracketing 
condition are shown in Figures 5.4 and 5.5. In Figure 5.4, the first two bracketing 
conditions are met, but F(B) < F(C) and so the minimum is not bracketed by the three 
points. Note that the minimum is in between B and C, but two points alone cannot satisfy 
the defined bracketing condition. In Figure 5.5, the full bracketing equation is met and the 
minimum is located within the interval [A,C]. 

The method that this routine uses is based upon parabolic curve fitting. Two initial 
guess points for A and B are supplied to the routine. If F(B) > F(A), then the values A and 
B are switched so that F(B) < F(A). The routine then calculates a guess for C such that the 
relative distance between A and B compared to the distance between A and C is the golden 
section, G r : 


C = B + (2 - G r ) (B - A) 


(5.26) 
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where: 


G r 


- 3 - V? 

2 


(5.27) 


After C is calculated, the routine does a parabolic curve fit between the three points. 
The routine solves for the minimum of this parabola and examines the function evaluated at 
this point. 

Examples of this procedure are shown in Figures 5.6 and 5.7 where points A, fi, and C do 
not bracket the function minimum. The abscissa of the parabolic minimum is U. In Figure 
5.6, the curve Fit results in U being located to the right of C. If the function to be 
minimized is F\, then the function value at U is less than than the function value at C. The 
minimum is not bracketed in this case. If Fi is the objective function, then the function 
value at U is greater than the function value at C. In this case, all of the bracketing 
conditions are met and the minimum is bracketed by B, C, and U. 

Figure 5.7 shows what would happen if the parabolic minimum was located between B 
and C. If the function is defined by Fi, then the function value at U is less than the 
function value at C. The minimum is bracketed by B, U, and C. If the function is defined 
by F 2 then the value of the function at U is less than B but still greater at C. The minimum 
of F 2 is not bracketed. 

If the function is not bracketed on the first iteration, then the procedure recalculates 
another U based on golden section magnification and does another parabolic curve fit with 
the three lowest points. The routine continues until the minimum has been bracketed. 

At this point, the method for choosing the initial points to give this routine will be 
discussed. Press proposed using the points A = 0 and fi = 1 as the initial test points. For 
the optimization of on-orbit mass, this initial guess for fi was frequently unsuitable. 
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As stated previously, the purpose of the one-dimensional minimization procedure is to 
find the value of t]k which when multiplied by the search direction, p*, and added onto the 
current estimate of the minimum, x*, would minimize the value of the objective function 
along the search direction at the new estimate, x*+i. The point A = 0 poses no problem 
because it translates into a test point for 77 * that will result in the current estimate, x* + i, 
being set to the previous estimate, x*. However, the point B = 1 corresponds to adding the 
entire search vector onto the current estimate. If the search vector is large, then the new 
estimate that would result from this addition could lie outside of the physical constraints of 
the trajectory. For example, the attitude at the end of launch maneuver cannot be larger 
than 90*. 

To correct for this problem, the following procedure was developed to find a suitable 
initial guess for 77 *. Suppose the constraints on the estimate x are specified such that each 
component of x has both an upper and a lower bound: 


li < Xi < Ui for i = 1 :n 


(5.28) 


where: 


/, = lower bound on x, 

Ui = upper bound on x, 

A value of 77 * is desired such that when it is used in equation 5.22, none of the 
components of the new estimate, x*+i, will violate the condition given by equation 5.28. 
To find such a 77 *, 2 n cases of equation 5.22 can be set up which would take one of the 
components of the current estimate to one of its constraints. For simplicity, the k subscript 
has been removed. A set of n equations can be defined to take the components of the 
estimate to their lower bounds: 


/, = xi + 77 pi for i = 1 :n 


(5.29) 


r O- 
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and another set of n equations can 
estimate to their upper bounds: 


be defined which would take the components of the 


Ui=Xi + T}pi for i = \:n 


Each equation will yield a unique value of TJ. There exist 2 n values of 77 : 


Vdesired € {T ]j) iOT j = 1: 2 H 


(5.30) 


(5.31) 


Of the 2 n possible values of 7 ] that could be found from the above equations, only n 
positive values exist . This is because the search direction will point towards either the 
upper or lower bound for each component of the estimate. Therefore none of the n 
negative values of rj can be chosen as the initial test point given to the bracketing routine 
because the optimal t\ used in equation 5.22 must be positive. 


T] desired e { T]j > 0} fmj = 1: 2/1 


(5.32) 


The value of t ] desired must finally equal the minimum of the above set divided by the 
quantity that the bracketing procedure will magnify it by to calculate the first guess of the 
third point of the bracketing interval, C: If A = 0, then the first guess of the quantity, C, is 
defined by equation 5.26 as: 


C = B + (2 - G r )B = (3 - G r )B 


(5.33) 


Since the initial guess for B could be magnified by (3 - G r ): 

minimum { rj; > 0 } 

Vdesired = 3^cT r f or J = 1: 2n (5.34) 

This value of T] guarantees that the initial triplet of abscissas used in the bracketing 
routine will not violate the constraints imposed on the estimate of the minimum. For the 
A.L.S. trajectory, the constraints on the trajectory parameters were chosen as: 
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75* < 6 f < 90' 


(5.35) 


T< ai < 17' 

(5.36) 

1*< a 2 < 17* 

(5.37) 


After the bracketing routine has found a triplet values of rj which satisfy the bracketing 
constraints, the set is input into a one-dimensional minimization routine. This routine is a 
variation of Brent's method. Brent's method uses a combination of the golden section 
search and the parabolic curve fitting procedure described previously for the bracketing 
routine. Brent's method, however, continues to narrow the bracketing interval until the 
minimum is found within a specified tolerance. The variation on Brent’s method used 
utilizes the derivative of the objective function with respect to T} k . 

An example of this is shown in Figure 5.8. The derivative of the function is evaluated 
at the middle point of the bracketing triplet, B. If the derivative is positive, the minimum 
will likely lie between A and B. In the case shown in the figure, the derivative is negative 
and the minimum of the function lies between B and C. The derivative at C is then 
evaluated. The two derivative values DF(B ) and DF(C) are linearly extrapolated to zero. 
The abscissa where the zero derivative occurs is used as the next trial point 

The method for calculating the derivative of the objective function with respect to rj ki 
DF(T] k ) is defined as follows: A new state estimate is calculated: 

x*+i = x* + Tj k p k (5.38) 

The gradient at this estimate is calculated: 

g* + i = VF(x* +1 ) (5.39) 
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Figure 5.8 Estimating Location of Function Minimum 
by Extrapolation (or Interpolation) of Slopes to Zero 


The derivative of the objective function with respect to 77 * is approximated as. 


DF{t)k) = 



- gLiP* 


(5.40) 
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5.5 Gradient Approximation 


The gradient vector of the objective function can be approximated using the "finite 
difference technique". For a sufficiently small scalar 8j, the Taylor Series can be used to 
approximate to first-order a function of n variables as: 

F(x + 6jtj) * F(x) + 5jgj(x ) (5.40) 


where: 

e ; = unit vector in the y-th coordinate direction 
gj = y-th component of the gradient vector 


The y'-th component of the gradient can then be approximated as: 


, x F{x + Sfij) - Fix) 

g. (x) = u- 


(5.41) 


This form is called the "forward difference approximation" and is exact only for a linear 
function. Another form can be derived by defining: 

Fix - Sjej) - Fix) - Sjgjix) (5.42) 


Subtracting of equation 5.42 from 5.40 yields: 


g;<x) - 


Fjx + Sfij) - Fjx - 6j€j) 

28 j 


(5.43) 


This form is called the "central difference" approximation" and is exact for a quadratic 
function. However, this form requires twice as many function evaluations as for the 
forward difference approximation. To reduce the number of computations involved in the 
optimization procedure, it was decided to approximate the gradient using the forward- 
difference approximation. 
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When computing the forward difference approximation of the gradient vector, the 
objective function, is first evaluated at the current state estimate x. Each component of the 
gradient vector is found by first perturbing the corresponding x component by a small 
amount, subtracting the function value at x, and dividing the difference by the perturbation. 
Thus, if an objective function is defined by n independent variables, then n + 1 function 
evaluations will have to be performed to approximate the gradient vector. 

The choice of the values to use for the perturbation 5 must take into consideration the 

accuracy of the predictive simulation and all of the subprograms it uses. If a value is 
chosen which is smaller than the precision of the simulation, then the value of on-orbit 
mass which is returned from the predictive simulation will not vary over the perturbation. 
Many sources, including Scales and Press, recommend using the square root of the 
machine precision as the perturbation value. This will only work if the procedure used to 
calculate the objective function has as high an accuracy as the machine. The perturbation 

used this study was chosen as: 

s=o.oor (5 ' 44) 

5.6 Conclusions 

This chapter compared various numerical optimization methods that could be used to 
automate A.L.S. trajectory design. The objective of trajectory design is the maximization 
of on-orbit mass. Because the optimization methods studied are designed to minimize a 
given objective function, the objective function is defined as the negative of the on-orbit 
mass . The conjugate gradient method was chosen over direct search methods because it 
utilizes first derivative information and thus has faster convergence properties. The 
conjugate gradient method was found to be more practical than Newton's methods for the 
on-orbit mass optimization problem because it did not require knowledge of the Hessian 
matrix. This matrix would be take too many function evaluations to approximate with finite 

differencing. 
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The chapter then outlined the conjugate gradient method and described the algorithm 
flow. The conjugate gradient method takes an initial guess for the state vector that will 
minimize the objective function and iterates on this guess using weighted search directions 
until the minimum has been located to within some specified tolerance. For a quadratic 
function of n variables, these search directions are mutually conjugate and will locate the 
minimum within n iterations. 

The one-dimensional minimization along each search direction was also discussed. 
Two routines that are given in Press were used. The first brackets the minimum with a 
triplet of points such that the function value at the central point is less than the function 
value at both of the other points. The second routine takes this bracketing set and continues 
to narrow the bracketing interval until the minimum is found to within some tolerance. 

A procedure was developed for making an initial guess for one of the bracketing points, 
i.e. the weight that is applied to the current search vector in the calculation of the new 
estimate. Too large a weight for the search vector would result in the updated estimate 
being located outside the physical constraints of the problem. 

Finally, the finite difference technique used to approximate the gradient vector was 
discussed. The forward difference approximation was chosen over the central difference 
approximation because it did not require as many function evaluations. For an objective 
function of n variables, the forward difference approximation requires n + 1 function 
evaluations. 
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Chapter Six 


SIMULATION AND EVALUATION 


6.1 Introduction 

The optimization scheme described in the preceding chapter is used in conjunction with 
the predictive simulation to find the set of trajectory shape parameters which will maximize 
on-orbit mass for a given Qa limit. The predictive simulation is initialized with the current 

vehicle state and is provided with a pre-launch wind velocity measurement. 

This chapter first describes in Section 6.2 a decision process used to define a Qa limit 
for pre-launch trajectory optimization. The chosen Qa limit will allow the vehicle to 
achieve a desired on-orbit mass in the presence of large in-flight wind dispersions from the 
measured pre-launch wind profile. 

Section 6.3 shows the ability of the optimization procedure to define optimal trajectories 
prior to launch for the trajectory shape described in Chapter 3. However, the use of this 
trajectory shape made the optimization process sensitive to the initial guess of the optimal 

trajectory shape parameters. 

This problem lead to the definition of a new trajectory shape which is described in 
Section 6.4. Pre-launch trajectory optimization results are presented for the new trajectory 
shape. The use of this trajectory shape results in an optimization process which is much 

less sensitive to the initial guess of the optimal solution. 

Section 6.5 discusses the problems encountered in attempting to do trajectory 
optimization over the flight shape parameters in flight. It was found that, using either 
trajectory shape, the only option for redesigning the trajectory in flight was to use a 
different value of the Qa limit than was used in pre-launch trajectory optimization. 
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6.2 Decision Process for Choice of Qa Limit 

The objective of this decision process is to establish a Qa limit to be used for trajectory 
design- both pre-launch and in flight. This Qa limit must be chosen so that the vehicle will 
meet a baseline on-orbit mass, mf in the presence of the largest wind dispersions 

from the pre-launch wind measurement that could occur in flight. The baseline on-orbit 
mass will be specified by mission planning and should include the dry mass of the core 
stage, the payload, the fuel mass needed for any maneuvers in orbit to deliver the payload, 
and enough fuel mass to accommodate any problems with the vehicle's insertion into the 
desired orbit. The vehicle will be able to accomplish this objective by carrying a fuel pad to 
account for the maximum wind dispersions. The amount of fuel pad needed is defined by 
mission planning. This decision process assumes that a bound exists for the maximum in- 
flight wind dispersions from the pre-launch wind measurement. 

The other assumptions made by this process are shown in Figure 6.1. This figure 
shows a hypothetical set of curves corresponding to the maximum on-orbit mass, m f , that 
could be obtained (using a trajectory optimization procedure) for a specified Qa limit for 
different wind profiles. The different wind profiles shown are for a strong headwind, 
HW max , a strong tailwind, TW max , and no winds, NoW. The amount of on-orbit mass 
that can be obtained will increase as the winds go from strong head winds to strong 
tailwinds. Headwinds oppose the vehicle's motion and degrade performance while 
tailwinds tend to improve performance. 

The loci of mf points are shown as curves that flatten out with increasing Qa limit 
because there is an upper bound on the on-orbit mass that can be achieved with increasing 
Qa limit. For a large Qa limit, the on-orbit mass will approach a constant value because 
the normal force constraint will, in effect, be removed. 
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Figure 6.1 : Assumptions Made by Decision Process 


The abscissas of the Qa limit axis shown in Figure 6.1 are: 

Qoctim [design spec ] = Qa limit specified by vehicle design 

The Qcc experienced by the vehicle in flight cannot exceed this 
value or normal loads on the vehicle could cause structural 
damage. 

Qaiim fn g ht\ - Q a lirnit to be used in flight for Qa limiting mode switch in Phase 3 

This limit must be less than the limit set by vehicle design by a 
certain factor of safety. 

Qaiim max ^ = maximum Qa limit to be used for trajectory design 

This limit must be less than the limit used for the in-flight Qa 
limiting mode switch to accommodate excursions in angle of attack 
caused by the estimator performance, noise, wind spikes, etc. 
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Q a iim[TD min] = minimum Qa limit to be used for trajectory design 

This limit sets a lower bond on the Qa limit that can be used for 
trajectory design. A smaller Qa limit than this value will not allow 
the vehicle to reach orbit. This value will increase in the direction 
from tailwinds to headwinds. 

Having explained the above concepts, it is now possible to outline the steps for the 
decision process: 

1. For the family of curves of on-orbit mass vs Qaum, one of these achieves an on- 
orbit mass equal to m/ baseline for a Qau m equal to Qaiun [jd max ]■ Th e headwind profile 

assumed for this curve is the maximum headwind that is acceptable to accomplish mission 
objectives. If the winds measured pre-launch are (100-S)% of this value or higher, where 
S% = percentage of HW max that bounds wind dispersions in flight, then the mission must 
be cancelled. 

2. The reader is now referred to Figure 6.2. If the pre-launch wind measurement is 
R% of HW max , then the maximum head winds that could be experienced in flight is 
(R+S)% of HW max . 

3. Locate the intersection of the (R+S)% curve and m / beeline horizontal line. The Qa 
limit where this occurs ( Qai mTD ) is the Qa limit to be used for trajectory design for an 

assumed wind of R% of HW ma x- 

4. The value of the on-orbit mass that can be achieved by the R% wind profile at 

QttlimrD m f TDinif 

5. The initial fuel pad to account for the maximum headwind wind disturbance, S% of 
HW max is then 


fPi M f JDinit ' m f baseline 


( 6 . 1 ) 
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NoW 



Figure 6.2: Decision Process for a Headwind Pre-Launch Measurement 

The fuel pad to be used during endoatmospheric flight should decrease linearly with 
time to be zero when an altitude is reached where the wind velocity is zero: 

< 6 - 2 > 

where tf Wiru i s — the approximate time when the wind velocity has decreased to zero. With 
this function for the in-flight fuel pad defined, the target mass of any trajectory design 
updating in flight should be: 

m f baseline + fPW ^ 
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63 Pre-Launch Trajectory Optimization 


Four sets of pre-launch trajectory optimization runs were carried out. They utilized the 
conjugate gradient method in conjunction with the predictive simulation run with a 1.0 
second integration time step. The first two sets corresponded to the angle of attack 
trajectory design profile described in Chapter 3. The second two sets correponded to a new 
angle of attack profile that was developed in the course of this study. The first two sets of 
runs will be referred to as the "old alpha profile" runs and the second two sets will be 
referred to as the "new alpha profile" runs. 

The first set of runs for the old alpha profile used 60% of the Vandenburg 69 (Van69) 
headwind profile and the second set used 100% of Van69. Tables 6.1 and 6.2 show the 
optimization results of these runs for different values of the Qa limit. The "major 
iterations" column refers to the number of different conjugate gradient search directions that 
were required for the on-orbit mass to converge to within a tolerance of 0.5 slugs. The 
"avg. ID iterations" column refers to the average number of iterations required for the one- 
dimensional optimization routine to find the scalar weighting of each search direction. For 
each run, the maximum on-orbit mass is given along with the values of the optimal set of 
trajectory shape parameters. Figure 6.3 shows a plot of the maximum on-orbit mass 
achieved on each run vs. the Qa limit. As expected, the maximum on-orbit mass increases 
for increasing Qa limit. In addition, the maximum on-orbit mass at each Qa limit is 
smaller for the 100% Van69 wind profile than for the 60% Van 69. Within the range of 
Qa limits tested, the plots are both approximately linear. 
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Qa Limit 

Major Iterations 

Avg. ID 
Iterations 

Maximum m/ 
(slugs) 

Optimum 
6 f , a\, a 2 (deg) 

2250 

3 

3 

11061.5 

IRERIHH 

2500 

7 

4 

11068.3 


2750 

4 

2 

11075.8 


3000 

8 

5 

11082.1 



Table 6.1: Old Alpha Profile Pre-Launch Optimization Results 
for 60% Van69 Headwinds 


Qa Limit 

Major Iterations 

Avg. ID 
Iterations 

Maximum m/ 
(slugs) 

Optimum 
e f , a u a 2 (deg) 

2500 

3 

5 

11052.9 


2750 

5 

3 

11059.5 


3000 

3 

4 

11066.8 



Table 6.2: Old Alpha Profile Pre-Launch Optimization Results 
for 100% Van69 Headwinds 
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Figure 6.3: On-Orbit Mass Plot for Old Alpha Profile 


Figure 6.4 shows the plots of angle of attack versus time for the optimal solutions of 
trajectory parameters for the case of 60% Van69 headwind. A similar set of plots is shown 
in Figure 6.5 for the 100% Van69 headwind case. This condition introduced an unforseen 
difficulty in the implementation of the optimization procedure. If the initial guess for the 
trajectory parameters was not chosen judicially, the procedure could drive the solution to a 
region in which the time spent in the a, region approached zero. At this point the gradient 
component corresponding to a, becomes identically zero and the assumptions underlying 
the optimization become invalid. It sometimes required that several different initial guesses 
for the parameters be considered before a solution was obtained in which the time spent in 
the a, region was greater than zero. 
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To counteract this problem, a new trajectory shape (alpha profile) was developed. The 
difference between the new alpha profile and the old alpha profile is in the launch maneuver 
(Phase 2) part of the trajectory design. Instead of steering the vehicle to a constant angle of 
attack at the end of the launch maneuver, the objective of the launch maneuver trajectory 
design was changed to steering the vehicle to the angle of attack limit that was defined by 
the Qa limit. The error for the one-dimensional launch maneuver optimization was simply 
based on Qa limit instead of the angle of attack. Table 6.3 shows the error in the Qa at the 
end of the launch maneuver. The error is a very small percentage of the Qa limit for each 
run. Instead of spending a very small amount of time in the constant alpha mode, the 
vehicle went directly to the angle of attack limit. The number of trajectory parameters was 
reduced from three to two because a\ was eliminated. This change in trajectory shape 
made it possible to input almost any realistic combination of 0/and ai into the 
optimization scheme and have the on-orbit mass converge to within the tolerance. In 
addition, the reduction in the number of trajectory parameters reduces the number of 
independent variables used in the optimization scheme. 


% Van69 Winds 

Qa Limit 
(lbs deg/ft**2) 

Error in Qa 
(lbs deg/ft* *2) 

60 

2250 

4.0 

60 

2500 

6.3 

60 

2750 

0.7 

60 

3000 

1.7 

100 

2500 

6.9 

100 

2750 

4.6 

100 

3000 

5.2 


Table 6.3: Error in Qa at End of Phase 2 for Different Simulations of 
the New Trajectory Design Method 
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The results for the new alpha profile are presented in the same way as for the old alpha 
profile for comparison between the two. Tables 6.4 and 6.5 show the optimization results 
for the new alpha profile for 60% Van69 headwinds and 100% Van69 headwinds. The 
maximum on-orbit mass obtained for each Qa limit is very close to those values obtained 
for the old alpha profile. This can be seen in Figure 6.6. The resulting alpha profiles for 
the optimal solutions are shown in Figures 6.7 and 6.8. 


Qa Limit 

Major Iterations 

Avg. ID 
Iterations 

Maximum m/ 
(slugs) 

Optimum 
Of, oil (deg) 

2250 

7 

5 

11060.7 

88.2, 7.8 

2500 

3 

6 

11067.5 

86.8, 9.6 

2750 

5 

6 

11075.0 

85.3, 11.3 

3000 

3 

4 

11082.0 

83.7, 9.6 


Table 6.4: New Alpha Profile Pre-Launch Optimization Results 
for 60% Van69 Headwinds 


Qa Limit 

Major Iterations 

Avg. ID 
Iterations 

Maximum mj 
(slugs) 

Optimum 
6 { , a 2 (deg) 

2500 

7 

5 

11052.7 

89.1, 10.1 

2750 

6 

3 

11060.3 

87.5, 10.0 

3000 

5 

4 

11066.6 

86.0, 9.2 


Table 6.5: New Alpha Profile Pre-Launch Optimization Results 
for 100% Van69 Headwinds 
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Alpha (deg) On-Orbit Mass (slugs) 



Q- Alpha Limit (lbs deg/ft**2) 

Figure 6.6: On-Orbit Mass Plot for New Alpha Profile 



Figure 6.7: New Alpha Profile: Alpha Plots for Optimal Solutions of 60% Van69 
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6.4 In Flight Trajectory Design 

It was decided that any trajectory design in flight should occur after the vehicle has 
completed the launch maneuver. During the launch maneuver, the altitude is still low 
enough such that the winds have not built up enough to cause the vehicle to deviate from its 
trajectory by a significant amount. 

Examining the prelaunch optimization results presented in the previous section, the 
optimal trajectories of both shapes spent most of their time flying the angle of attack limit in 
the constant Qa phase. Thus, in-flight trajectory optimization over a x and a 2 is not 
practical because the trajectory is almost independent of both parameters. In addition, as 
the predictive simulation is begun later in flight, even less time is spent within the a 2 

portion of flight. This trend was illustrated in Chapter 4. 

Therefore, the only redesign that can occur in flight with the trajectory shapes as 
defined, is based on varying the Qa limit. The prelaunch optimization results proved that 
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the on-orbit mass would increase with increasing Q a limit. This can be used to ensure that 
the vehicle meets its time-varying fuel pad requirement, fp(t). 

As an example of this, a run was made where the winds experienced in flight were 
headwinds stronger than those measured before launch. Referring again to Figure 6.3, a 
baseline on-orbit mass was specified as 11060 slugs. A wind measurement taken before 
launch was assumed to be the 60% Van69 wind profile. The winds were assumed to be 
able to vary as much as 40% of Van69. The Qa limit where the 100% Van69 curve 
crosses the baseline mass was 2750 lbs/ft 2 . The on-orbit mass of the 60% Van69 curve at 
this Qa limit is 1 1076 slugs. The initial fuel pad is therefore defined as 16 slugs. 

The set of optimal trajectory parameters for 60% Van69 at a Qa limit of 2750 lb-deg/ft 2 
was then used to design a trajectory for the full simulation. The full simulation was run 
using all of the estimators. The winds used in flight were 100% Van69. Without any 
trajectory redesign in flight, the on-orbit mass of the vehicle was 1 1052 slugs. 

Another full simulation was run under the same conditions, but at 60 seconds the 
trajectory was redesigned. At 60 seconds, the predictive simulation was called with the 
current vehicle state, the prelaunch wind profile, and a Qa limit of 2750. The predictive 

simulation returned an on-orbit mass of 1 1063 slugs. At 60 seconds, the fuel pad will have 
decreased linearly to 8 slugs. The baseline on-orbit mass plus this fuel pad equals 11068 
slugs. Therefore, a search was made to find the Qa limit that would meet this requirement. 
A Qa limit of 2875 was found to satisfy the requirements and the trajectory was redesigned 
using this Qa limit. The vehicle flew from 60 seconds on with the new trajectory. The on- 
orbit mass resulting from this change was 11055 slugs. The performance of the full 
simulation with and without trajectory design is shown in Figure 6.9. 

An engine out case was attempted, but the vehicle was not able to achieve the desired 
orbit, even with trajectory design in flight at a higher Qa limit. 
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Chapter Seven 


CONCLUSIONS AND RECOMMENDATIONS 

7.1 Conclusions 

One of the objectives of this thesis was to develop a computer program that would 
automate the prelaunch trajectory design process of the Advanced Launch System 
(A.L.S.). The other objective was to explore the possibility of redesigning the trajectory in 
flight to compensate for wind dispersions or an engine out failure. 

Since the objective of the A.L.S. trajectory design was the maximization of on-orbit 
mass, a numerical optimization scheme was employed to determine the set of trajectory 
parameters, for a given trajectory shape, which maximize on-orbit mass. A predictive 
simulation was developed which could accurately calculate the on-orbit mass given the 
vehicle state, wind information, and a set of trajectory parameters. This simulation utilized 
an idealized control assumption that the control variable, either pitch attitude or angle of 
attack, was perfectly achieved within a given integration time step. Because of this 
assumption, the rotational equation of motion did not have to be solved to find pitch attitude 
and a larger integration time step could be used. Thus, the computational speed of the 
predictive simulation over the full simulation was shown to be close to that of the full 
simulation, even with integration time steps as large as one second. 

After reviewing the current literature, a conjugate gradient method was chosen over 
other types of numerical optimization schemes. The conjugate gradient method requires 
only knowledge of the objective function and its gradient. The objective function was 
defined as the negative of on-orbit mass. The gradient is the vector of derivatives of the 
on-orbit mass function, calculated using the predictive simulation, with respect to the 
trajectory parameters. The gradient is approximated with finite differencing. 
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The optimization procedure was applied to the problem of prelaunch trajectory design. 
Using the trajectory shape as originally proposed, it was discovered that optimal solutions 
produced trajectories in which the time spent in the constant angle of attack portions of 
Phase 3 was very small. As a consequence, it was found that optimal solutions were very 
dependent on the initial guess for the trajectory parameters. In some cases, when the initial 
guess was not close to the optimal solution then convergence to a realistic solution could 
not be obtained. In these cases, the solution was driven to a condition in which the time 
spent in the first constant angle of attack portion of Phase 3 was identically zero. In this 
case, the assumptions underlying the optimization process became invalid and a correct 
solution was not obtained. 

To correct this problem, a new trajectory shape was defined which was similar to the 
old shape except that instead of having the launch maneuver take the vehicle to a constant 
angle of attack at the start of Phase 3, the vehicle was taken directly to the angle of attack 
defined by the Qa limit. With this trajectory shape, neither component of the gradient ever 
became identically zero before convergence was obtained. The procedure was therefore 
able to converge to the optimal solution. In addition, the change of trajectory shape 
reduced the size of the parameter set from three to two and, consequently reduced the 
dimensions of the optimization problem. Even with the modification in trajectory shape, 
the predicted on-orbit masses for the optimal solutions were approximately the same as 
those obtained using the old trajectory shape. 

Both sets of prelaunch optimization results from the old and the new trajectory shapes 
showed the vehicle spending a small amount of time in the second constant angle of attack 
portion of Phase 3. As a consequence, it was impractical to reoptimize for this parameter 
during flight. Since the optimum solution tended to approach a simple constant Qa 
trajectory after the launch phase, the only practical redesign in flight was to use the 
predictive simulation to design a new trajectory based on a different Qa limit. If it was 
determined from the predictive simulation that the on-orbit mass would be lower than that 
needed to meet the fuel pad requirement at the redesign time, then the trajectory would be 
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redesigned with a higher Qa limit. If it was found from the predictive simulation that the 
mass would be higher than needed to meet the fuel pad requirements, then the trajectory 
could be redesigned with a lower Qa limit to reduce stress on the vehicle. 

7.2 Recommendations 

Further study should be made to determine if more effective trajectory optimization 
parameters could be substituted for the constant angle of attack parameters of Phase 3. In 
particular, the constant angle of attack parameter at the end of Phase 3 might be replaced by 
a parameter that provides a more efficient transition to the Powered Explicit Guidance 
Phase (PEG). The choice of a new parameter for this transition phase could improve the 
effectiveness of the in-flight trajectory optimization. 

The launch maneuver design for this thesis uses a specified pitch rate profile of 
sinusoidal form which the vehicle is commanded to follow. Other launch maneuver 
profiles should be investigated that might provide improvements in payload capabilities. 

The impact of an engine-out failure was considered in this thesis, but it was found that 
it was not feasible to achieve the desired on-orbit mass for the assumed jet model. It was 
assumed, however, in this thesis that no throttle-up capability existed for the non-failed 
engines. Alternative designs should be investigated in which the engines nominally thrust 
at less than full thrust and then increase thrust in the event of an engine-out. For this 
assumption, the on-orbit mass requirements could more easily be met. It should be 
determined to what extent the trajectory and on-orbit mass are affected by the particular jet 
that is failed. A failed jet with a large moment arm from the longitudinal axis through the 
eg might require a cant angle bias of the remaining jets to reduce the average moment. This 
could in turn tend to produce a different angle of attack profile than in the nominal case. 
This occurrence might require modification of the trajectory optimization algorithm. 
Control problems associated with an engine-out should also be investigated. In particular, 
the ability to respond to transients produced by the engine-out should be studied. 
Provision should be made to modify control gains to maintain stability margins for the new 

thrust model. 
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Another optimization procedure could also be developed to automatically re-tune the 
gains used in the vehicle's control systems given last minute changes in payload, vehicle 
configuration, or mission objective. The objective function could be described in terms of 
the deviation from desired stability margins. This automation would save much in the way 
of prelaunch mission preparation. 

For this thesis, the steering method employed for Phase 3 was acceleration direction 
steering. A study should be made to compare the effectiveness of the trajectory 
optimization procedure with other steering methods e.g., flight path angle steering. 
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appendix 


/ 

NASA Langley Research Center provided this study with wind measurements taken at 
Vandenburg AFB, CA. The wind information was given as velocity plotted versus 
altitude. The profile used in this thesis is #69, shown in Figure A.l. For simplification, 
the profile as linearized as shown in Figure A. 2. 



Figure A.l Vandenberg #69 an #70 wind profiles. 
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Figure A. 2 Linearized Vandenberg #69 and #70 wind profiles. 
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